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SECTION  I 
INTRODUCTION 


1.1  THE  PROBLEM 

Within  many  subtle  limitations,  laser  velocimeter 
(LV)  systems  provide  non- interfering  point  measurements  of 
a  velocity  field  when  suspended  scatterers  pass  through  a 
small  optical  probe  volume  [1] .  For  wind  tunnel  measurements 
with  little  or  no  particulate  seeding  of  the  gas,  the  LV 
systems  of  concern  here  make  essentially  instantaneous 
sample  measurements  from  sparsely  distributed  scatterers 
which  are  detected  at  random  intervals  in  time.  The  primary 
objective  of  the  research  reported  herein  has  been  the 
development  of  digital  statistical  data  processing  tech¬ 
niques  for  the  estimation  of  turbulence  intensity,  auto¬ 
correlation,  and  power  spectrum,  with  low  mean  data  rates 
which  minimize  the  particulate  seeding  requirement.  An 
additional  objective  has  been  the  reduction  of  system  costs 
by  elimination  of  the  need  for  expensive  mass  memory  through 
the  use  of  real-time  minicomputer  processing. 

1.2  BACKGROUND  FOR  STATISTICAL  ESTIMATION 
1.2.1  Conventional  Spectrum  Analysis 

Many  engineers  who  are  familiar  with  the  use  of 
analog  spectrum  analyzers  in  measuring  periodic  data  are 
unaware  of  the  many  subtleties  in  the  measurement  of  power 


1 


AEDC-TR -74-53 


spectra  of  continuous  noise  processes.  However,  theory 
and  practices  are  well  developed,  at  least  for  stationary 
Gaussian  processes  [2,3]. 

In  many  applications  the  need  for  greater  accuracy 
or  the  digital  nature  of  the  original  data  format  has  led 
to  digital  power  spectrum  analysis  of  periodically  sampled 
data.  There  is  a  clear  distinction  between  the  require¬ 
ments  of  periodic  or  nearly  periodic  signals  and  those  for 
broadband  random  signals,  but  the  techniques  and  limitations 
are  well  understood  for  stationary  Gaussian  signals  [2-9] . 

An  excellent  introduction  to  this  subject  has  been  provided 
in  the  very  readable  paper  by  Richards  [4] .  Free  use  of  this 
background  material  is  utilized  in  this  report. 

1.2.2  Randomly  Sampled  Signals 

Fourteen  years  of  rigorous  mathematical  research 
with  statistical  averages  or  infinite  time  averages  indi¬ 
cates  that  not  only  is  it  possible  to  make  mean,  correlation, 
and  spectral  estimates  using  randomly  timed  samples  but  also 
that  there  is  a  significant  advantage  in  doing  so.  For  esti¬ 
mating  second  order  statistical  parameters  such  as  correlation 
and  spectrum  it  is  well  known  that  with  periodic  sampling  the 
sample  rate  must  exceed  twice  the  signal  bandwidth  to  avoid 
aliasing  error  [5],  This  Nyquist  criterion,  which  is  also 
required  for  time  history  reconstruction  from  periodic  samples 
often  imposes  difficult  constraints  on  the  data  processing 
when  spectral  resolution  requirements  are  added.  It  has  been 
shown  mathematically  that  time  history  reconstruction  from 
random  samples  also  requires  the  Nyquist  criterion;  but  practi 
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cal  reconstruction  with  low  error  requires  approximately 
three  times  the  Nyquist  rate  even  for  nonrealizable ,  non¬ 
linear,  straight-line  interpolation  and  much  higher  rates 
for  optimum  linear  filters  or  sample-and-hold  interpolators 
[10,11].  In  contrast,  alias-free  spectrum  estimation  may 
theoretically  be  obtained  with  any  mean  sample  rate,  however 
small,  when  the  random  sampling  process  and  the  signal  are 
both  stationary  and  independent  [12-14]. 

The  economical  implications  of  the  mathematical 
theory  are  very  significant  for  many  applications  in  addition 
to  the  LV  problem  which  naturally  provides  randomly  sampled 
data:  we  need  only  to  reduce  the  mean  sample  rate  to  reduce 
the  required  electronic  computation  speed  and  memory  to  ob¬ 
tain  economical  real-time  processing. 

Unfortunately,  in  contrast  to  the  rigorous  mathe¬ 
matical  theory  based  on  infinite  averages,  the  theory  of 
practical  spectral  estimation  from  finite  data  sets  is  in 
its  infancy  and  is  almost  non-existant .  First  we  must 
distinguish  between  a  few  different  classes  of  problems. 
Singleton  [15]  has  considered  estimation  of  parameters  of 
randomly  samples  signals  known  to  be  periodic  and  Shaw  [16] 
has  considered  periodic  sampling  with  random  time  jitter; 
but  these  are  different  problems.  In  particular,  our  experi¬ 
ments  have  shown  that  randomly  sampled  periodic  data  imposes 
minor  requirements  compared  with  randomly  sampled  broad¬ 
band  noise  processes.  Jones  [17]  has  investigated  another 
class  of  problems,  which  is  more  closely  related  to  our 
concern,  where  independent  random  omissions  of  periodic  sam¬ 
ples  occur.  This  type  of  sampling  becomes  similar  to  Poisson 
random  sampling  as  the  periodic  sampling  frequency  increases 
and  the  probability  of  a  given  sample  occurring  goes  to  zero 
[11]. 
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1.2.3  Estimation  of  Broadband  Spectra  from  Random  Samples 

In  the  literature  search  only  a  few  recent  papers 
were  found  which  relate  directly  to  the  problem.  Adegbola 
[18]  has  theoretically  analyzed  one  approach  in  detail  in 
which  the  data  are  used  to  form  a  rectangular  approximation 
to  the  Fourier  integral  formulas  applicable  to  continuous 
data.  He  presents  no  experimental  work  and  the  estimator 
equation  is  not  amenable  to  rapid  computation.  We  have 
derived  a  related  estimator  which  involves  less  approxima¬ 
tion  and  computation.  This  estimator  and  an  example  of  the 
poor  experimental  results  which  were  obtained  with  it  for 
subNyquist  mean  sampling  rate  are  presented  in  Appendix  I 
as  a  warning  to  others.  We  did  not  investigate  Adegbola* s 
estimator  experimentally  because  it  appears  slow  computa¬ 
tionally. 

In  a  very  recent  research  study  concerning  spectral 
analysis  from  single-particle  LV  data,  Asher,  Scott,  and 
Wang  [19]  attempted  a  straightforward,  economical  analog 
approach  in  which  continuous  sample  and  hold  signals  were 
analyzed  by- conventional  spectral  analysis  instrumentation. 
Their  unfortunate  results  should  be  considered  by  those 
who  would  apply  the  output  of  a  frequency- tracker  LV  system 
which  incorporates  hold  circuits  during  signal  drop  out  to 
conventional  instrumentation.  It  is  not  difficult  to  see 
by  hindsight  that  the  use  of  sample-and-hold  is  mathematic¬ 
ally  equivalent  to  our  rectangular  numerical  approximation 
of  the  Fourier  integrals  previously  mentioned.  For  practical 
single-particle  LV  data  rates,  or  for  trackers  with  appre¬ 
ciable  drop  out,  we  end  up  spectrally  analyzing  a  random  step- 
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function  whose  low-pass  spectrum  obliterates  the  desired 
spectrum. 

Another  digital  approach  was  taken  by  Thompson 
[20]  whose  derivation  was  inadequate  to  predict  the  addi¬ 
tive  white  spectral  level  he  obtained  in  simulations,  and 
whose  simulations  with  deteriministic  data  are  misleading 
in  relation  to  the  actual  difficulties  encountered  with 
signals  which  are  broadband  random  processes.  Through  a 
different  and  more  exact  reasoning  process  we  have  derived 
an  estimator  similar  to  Thompson's  which  correctly  sub¬ 
tracts  the  additive  white  level.  Our  estimator,  which 
is  based  on  the  known  relationship  between  the  signal  spec¬ 
trum  and  the  spectrum  of  the  samples  signal  [11,  21]  and 
its  experimental  demonstration  are  also  displayed  in  Appen¬ 
dix  I.  This  estimator  was  not  pursued  because  it  appears 
less  amenable  to  real-time  computation  than  the  approach 
taken . 

The  spectrum  estimator  which  we  derived  and  found 
most  useful  during  the  exploratory  phase  of  this  research 
consists  of  determining  an  autocovariance  estimate  from 
the  data  after  discretizing  the  occurrence  time  intervals 
as  though  the  sampling  were  periodic  with  most  of  the  samples 
randomly  missing.  Our  conclusion  that  this  approach  is 
optimum  at  the  present  time  is  supported  by  additional  ref¬ 
erences  by  Jones  [22,23]  in  which  he  describes  three  estima¬ 
tors  including  the  discretizing  autocovariance  approach  and 
concludes  that  this  approach  is  the  most  practical  of  the 
three.  Jones  unfortunately  makes  the  erroneous  conclusion 
that  once  the  time  discretizing  has  been  accomplished,  the 
error  analysis  problem  reduces  to  that  of  randomly  missed 
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samples  and  references  his  earlier  analysis  [17] .  This 
is  unfortunate  because  the  analysis  referred  to  does  not 
correctly  predict  the  spectrum  estimate  error,  as  we 
have  shown  by  experiment. 

Further  support  of  the  discretizing  autocovar¬ 
iance  approach  may  be  found  in  experimental  development 
results  currently  being  obtained  at  Lockheed-Georgia  Co. 
under  Contract  No.  F-33615-73-C-2032.  At  the  time  of  this 
writing,  good  agreement  between  hot  wire  anemometry  and  LV 
spectrum  measurements  has  been  obtained,  but  the  results 
are  not  yet  published.* 

1 . 3  SCOPE 


In  Section  2,  we  present  equations  for  estimation 
and  error  prediction  of  turbulence  parameters  measurable  by  a 
single-point,  single -component  LV  system.  Section  3  pre¬ 
sents  experimental  evidence  which  supports  the  error  predic¬ 
tion  equations  of  Section  2.  Section  4  provides  a  guideline 
procedure  for  experiment  and  implementation  design  with  fur¬ 
ther  interpretation  of  the  error  prediction  equations.  Sec¬ 
tion  5  is  a  brief  discussion  of  the  results,  and  Section  6 
provides  conclusions  and  recommendations  for  further  work. 


*  The  results  have  been  described  in  presentations 
by  D.  M.  Meadows  at  the  ’’Workshop  on  Laser  Doppler  Velocimetry , " 
held  at  Purdue  University,  West  Lafayette,  Ina.,  March  '27-29, 
1974.  Papers  will  appear  in  the  Conference  Proceedings. 
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SECTION  II 

THEORY 


2.1  INTRODUCTION 

2.1.1  Assumptions 

For  purposes  of  analysis,  we  consider  an  idealized 
model  of  the  fluid  flow,  the  LV  system,  and  the  mathemati¬ 
cal  computations.  The  practical  significance  of  these 
assumptions  is  discussed  later. 

The  velocity  component  U(t)  is  an  ergodic  Gaussian 
random  process  with  mean1 


UttT  =  U 


(2.1) 


and  variance  or  mean- square  turbulence  intensity. 


u  (t)  =  (U(t)  -  U)2  -  a2 

The  autocorrelation  of  U(t)  is  a  real,  even  function 


R(t)  =  U(t)  U ( t+T )  =  C(t)  +  U2 


(2.3) 


where  the  autocovariance  C(t)  is 

C (t ) -.  =  u(t)  u(t+r)  (2.4) 

The  two  sided  spectrum,  S(f),  and  the  autocovariance  form 
a  Fourier  transform  pair 

s(f)  =  C(T)  exp  [ - j 2irfT ]  dr  (2.5) 


1The  horizontal  bar  denotes  statistical  expectation 
or  infinite  time  average  interchangeably. 
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The  LV  system  is  assumed  to  detect  error -free 
instantaneous  samples  of  the  velocity  at  time  instants 
t^  which  are  the  realization  of  a  stationary  Poisson  point 
process,  statistically  independent  of  U(t)  ,  with  mean  rate  X. 

2.1.2  Mathematical  Background 


For  the  Poisson  process  the  expected  number  of 
samples  occurring  in  any  interval  of  length  T  is  XT,  with 
variance  XT,  and  the  probability  of  m  samples  in  T  is 


P  (m,T) 


(XT)m 


(2.6) 


It  proves  convenient  at  times  to  divide  the  time 
axis  into  uniform  increments  of  length  At  which  are  small 
enough  that  XAt<<1.  In  this  case,  expansion  of  (2.6)  as  a 
power  series  shows  that  the  probability  of  samples  occurring 
in  any  At  interval  is 

P(0,At)  -  1-p 

P (1 ,  At)  =  p  =  XAt  (2.7) 

P(m>l,  At)  =0 


For  sufficiently  small  At,  a  negligibly  small  timing  error 
occurs  in  replacing  the  exact  occurrence  time  t^  with  jAr, 
the  location  of  the  nearest  uniform  grid  element.  Using 
this  approximation,  the  data  sequence,  discretized  in  time, 
is  given  by  U(j)  I(j)  where  U(j)  is  the  value  of  U(t)  at 
jAr,  and  I(j)  is  an  independent  binary  random  variable  with 


P{I=1}  =  p 
P{I=0}  =  1-p 


(2.8) 
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The  use  of  this  notation  is  only  justified  when  p  =  A At 
is  small. 

The  concepts  of  equivalent  power  bandwidth  B 
and  turbulence  time  scale  tc  will  prove  useful: 


2B  = 


C  S(f)  df 

S(fpeak> 


.  o' 


s  ^peak> 


(2.9) 


T 


C 


C(t)  dt 


(2.10) 


where  s(£peak)  is  the  maximum  value  of  the  power  spectrum. 
For  a  typical  low-pass  spectrum  with  maximum  value  at  zero 
frequency,  we  note  that  C(t)  is  even  and 


SCO) 


U 


C(T)dT 


(2.11) 


to  obtain 

■  1/4  (2.12) 

2.2  MEAN  AND  INTENSITY  ESTIMATORS 

2.2.1  Introduction 

Well-known  statistical  formulas  exist  for  esti¬ 
mation  of  mean  and  variance,  and  for  prediction  of  the  mean- 
square  error  of  these  estimates,  when  uncorrelated  samples 
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are  used  as  data.  However,  it  is  contradictory  to  assume 
uncorrelated  samples  obtained  in  a  finite  time  with  random 
times  of  occurrence  when  we  require  sufficiently  high  data 
rates  to  estimate  correlation  functions  I 

The  detrimental  effects  of  correlated  samples  are 
as  follows:  the  standard  unbiased  estimate  of  variance 
becomes  biased,  and  the  variance  of  the  estimates  increases 
due  to  the  reduced  information  content  of  the  data.  We  have 
obtained  approximate  formulas  for  the  estimate  errors  which 
include  finite-time  correlation  effects.  These  approximate 
formulas  were  obtained  by  assuming  that  in  the  limiting 
case  of  high  A  the  error  could  be  no  less  than  that  which 
would  be  obtained  from  a  continuous  time  segment  of  U(t). 

The  errors  for  continuous  estimators  were  derived  from  results 
in  Papoulis  [21]  and  added  to  the  standard  formulas.  The 
details  of  the  derivations  are  omitted  for  brevity. 

2.2.2  Mean  Flow 


Subject  to  our  initial  assumptions  of  section 
2.1.1  the  usual  sample  mean 


-  1  N 

U  =  ry  L  U  • 

N  i=i  i 


(2.13) 


is  unbiased  and  has  rms  error  given  approximately  by 


rr*  ~  /1T4AT-^ 

°u  ”  a  c) 


.  2-  (1+^)% 


( 2 .14) 


where  the  X/2B  term  results  from  finite  time  effects. 
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2.2.3  Mean- Square  Intensity 

Subject  to  the  assumptions  of  section  2.1.1,  the 
conventional  variance  estimator  may  be  used  if  limitations 
due  to  finite-time  effects  are  considered: 

*2  1_  N 

0  '  N-l  ?  ,  (UrU)!  (2.15) 

1*1 


With  N^XT  we  find  the  continuous  estimate  of 
variance  over  the  same  duration  would  be  biased  and  assume 
we  can  expect  no  better: 

a2  -  a2  (1  -  — -  a2  (1  -  (2.16) 

The  bias  error  should  generally  not  be  significant  but  could 
easily  become  so  under  some  of  the  implementation  schemes 
discussed  later  in  which  data  is  broken  into  small  segments 
for  processing.  In  such  an  implementation,  care  should  be 
taken. 

The  derivation  of  the  added  finite-time  variance 
of  the  mean- square  intensity  estimate  is  simplified  by 
assuming  the  bias  is  negligible.  When  combined  with  the 
standard  result  given  by  Jenkins  and  Watts  [9],  the  result  is 

(°2  '  -  RTT  (2.17) 

where 

1 

l'c  =  O*  f0  C2(T)dT<Tc  (2.18) 
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The  inequality  is  valid  for  mono tonic ally  decreasing  posi¬ 
tive  functions  C(x)  but  should  be  checked  in  more  general 
cases . 


2.2.4  Normalized  Turbulence  Intensity 


We  may  estimate  I  =  a/ U  using  the  previous 
formulas  as 


0 


(2.19) 


Assuming  the  errors  in  a2  and  U  are  small,  Taylor  series 
expansion  shows 


I  “  1(1  ” 


(2.20) 


Assuming  that  the  bias  error  is  negligible,  and  that  the 

A  A 

errors  in  a2  and  U  are  uncorrelated,  we  find  the  variance 
approximately  as 


d  -  u*  *  (tt  +  +  ^ 


(2.21) 


In  low  turbulence  measurements,  the  I4  term  will  be  negli¬ 
gible  in  comparison  with  the  I2  term. 
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2 . 3  AUTOCOVARIANCE 

2.3.1  The  Estimators 


The  method  we  have  found  most  useful  for  spectrum 
analysis  begins  with  a  discretized  autocovariance  esti¬ 
mate  similar  to  that  reported  by  Jones  [22,  23].  The  exact 
time  differences  between  pairs  of  samples  are  represented 
approximately  as  integral  binary  numbers  with  an  arbitrary 
time  increment  At.  For  each  pair  of  samples  for  which  the 
discretized  time  difference  does  not  exceed  the  preselected 
maximum  lag  number  M,  the  product  of  the  two  samples  (lag 
product)  is  accumulated  in  a  lxM  matrix  SUM(k)  and  1  is 
added  to  a  corresponding  histogram  H(k).  For  purposes  of 
analysis  we  assume  that  the  true  mean  has  been  previously 
subtracted  from  the  data.  After  N  samples  are  processed 
the  resulting  estimate  C(k)  of  the  value  of  C(t)  at  T=kAT, 
k^O,  is 

COO  ■  HTtT  SUM00  (2.22) 


which  the  summation  includes  only  the  H(k)  products  for 


-  t. 


At 


k  |  <  0.5 


(2.23) 


For  k=0,  the  practical  advisability  of  including  products 
for  which  |t^  -  tj j  <  At/2  along  with  the  u?  terms  depends 
on  the  implementation;  without  these  terms,  H(0)  =  N  and 


C(0) 


1=1 


(2.24) 
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A  very  similar  autocovariance  estimate  may  also 
be  expressed  in  the  alternate  notation  which  assumes  assign¬ 
ment  of  t^  to  the  nearest  value  of  jAx  prior  to  computation 
as 

N'-k 

l  u(j)  u(j+k)  I(j)  I(j+k) 
j=0 

C'(k)  »  - —  (2.25) 

N'-k 

I  Kj)  I(j+k) 
j  =  0 


when  N=sAT=AN'At  and  I(j)  is  the  binary  indicator  sequence. 
The  estimator  is  denoted  with  a  prime  to  remind  us  that 
(2.25)  differs  from  (2.23)  and  (2.24)  because  ±At/2  errors 
in  the  discretization  of  t^  and  t^  prior  to  subtraction 
allows  maximum  errors  in  time  differences  of  ±Ax.  The  max- 

a  /v 

imum  time  error  for  C(k)  is  ±Ax/2.  Also  C(k)  uses  a  fixed 
number  of  samples,  while  C'(k)  uses  a  fixed  length  of  time. 
The  relation  N  =  AT  is  only  true  on  the  average. 

2.3.2  Bias  Error  Due  to  Time  Increment 


Jones  dismisses  the  analysis  of  the  discretized 
autocovariance  estimate  by  reference  to  his  analysis  of 
periodic  sampling  with  random  omissions  [17],  If,  in  fact, 
Ax  is  assumed  small  enough  to  make  the  time  errors  negli- 
gible  in  both  C(k)  and  C'(k)  then  the  application  of  condi¬ 
tional  expectation  techniques  to  the  random  omission  model 
shows  that  the  estimates  are  unbiased  provided  that  N  is 
sufficiently  large  so  that  each  lag  value  has  at  least  one 
lag  product.  We  note,  however,  that  practical  limitations 
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will  sometimes  require  Ax  to  be  as  large  as  is  permissible 
without  bias  error.  In  the  case  of  periodic  sampling  with 
random  omissions  Ax  must  satisfy  the  Nyquist  criterion. 

In  Appendix  II  we  demonstrate  that  the  expected 

A 

value  of  C(k)  is 


CCk)  *  /  C(x)  pk(x)dx  ( 

where  pk(x)  is  the  conditional  probability  density  for 
occurrence  of  lag  delays  and  is  nonzero  only  on  the  interval 
(kAx  -  Ax/2,  kAx  +  Ax/2).  The  function  pk(x)  has  unit  area 
and  is  not  given  by  a  simple  formula.  The  somewhat  compli¬ 
cated  exact  expression  for  pk(x)  is  given  in  Appendix  II 
for  the  case  of  uniform  Poisson  random  sampling. 

As  we  show  in  section  2.3.4  below,  in  many  prac¬ 
tical  cases  pk(x)  may  be  replaced  by  a  constant  1/Ax  over 
the  interval  of  definition  with  the  result  that  the  expected 

/v 

value  of  C(k)  is  the  average  value  of  C(x)  on  the  surround¬ 
ing  ±Ax/2  interval. 

The  limitation  imposed  by  selection  of  Ax  is 
best  understood  by  its  effect  on  the  spectrum  estimate  and 
will  be  discussed  further  in  section  2.4.  The  resulting 
maximum  value  of  Ax  is  of  the  same  order  of  magnitude  as 
the  Nyquist  criterion  for  periodic  sampling.  The  bias 

A  A 

error  for  C'(k)  is  expected  to  be  worse  than  that  for  C(k). , 
but  on  the  same  order  of  magnitude. 


.26 
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2.3.3  Additional  Bias  Error 


For  the  present  we  assume  the  true  mean  U  has 
been  subtracted  from  the  samples  before  formation  of  the 

A 

estimate  C(k).  In  implementation  it  is  actually  necessary 
to  subtract  sample  means  to  avoid  low-frequency  distortions 
of  the  spectrum  estimate.  This  procedure  is  expected  to 
produce  small  bias  errors  in  the  autocovariance  estimate 
which  are  associated  with  removal  of  the  very  low-frequency 
components,  but  the  effect  has  not  yet  been  analyzed  in 
detail . 


2.3.4  The  Histogram 


The  histogram  H(k)  of  the  number  of  lag  products 

obtained  at  each  lag  value  k  plays  a  significant  role  in 

/\ 

practical  application  of  the  estimator  C(k).  In  Appendix 
II  we  have  evaluated  its  expectation  approximately  under 
the  assumption  that  XAt  is  small  and  N  is  large  with  the 

A 

result  that  with  only  u?  terms  used  in  C(0), 


H(k=0)  =  N 


(2.27) 


H(k)  *  NXAt  -  k(XAx) 2  ■  (^1  -  k)  (XAt)2 


The  result  shows  that  the  number  of  products  at  each  lag 
is  generally  much  less  than  the  number  of  u?  terms.  The 
expected  number  of  products  at  each  lag  decreases  linearly 
with  k,  but  if  the  data  collection  time,  T,  for  each  data 
segment,  is  much  greater  than  the  maximum  delay,  i.e.,  if 

T  >>  MAt,  N  >>MXAt  (2.28) 
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then  the  value  of  H(k)  is  nearly  constant.  This  condition 
is  an  important  consideration  in  implementation  which  must 
be  observed  when  segmented  data  methods  are  used. 

The  variance  of  H(k)  has  not  been  theoretically 
evaluated,  but  we  expect  it  to  behave  in  the  manner  of  a 
Poisson  process  for  acceptably  large  values;  i.e. 

(HCk)  -  H(k))2  -  H(k)  f2.291 


We  may  obtain  an  approximate  expression  for  the 
function  pk(x)  utilized  in  section  2.3.2  from  FT(k)  by 
assuming  that  H(k)  is  an  unnormalized  estimate  of  the 
probability  density  for  occurrences  of  lag  values  on  the 
range  0<t<MAt.  The  result  is 


p  (t)  *  *** 

NXAt  -  kX*Ax2 


(2.30) 


It  is  expected,  therefore,  that  pk(x)  may  be  replaced  in 
(2.26)  by  a  unit  area  rectangular  function  of  width  Ax 
and  height  1/Ax  when  the  data  collection  time  is  much 
greater  than  the  maximum  delay. 

2.3.5  Variance  of  the  Autocovariance  Estimate 


When  the  system  parameters  are  correctly  chosen 
to  avoid  bias  error,  the  important  remaining  problem  is 
the  determination  of  the  statistical  error  of  the  estimate 
in  terms  of  the  amount  of  data  collected  and  processed. 
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In  Appendix  II  we  derive  the  variance  of  C(k)  under  the 
assumptions  that  the  normalized  rate  parameter  X/2B  is 
small  compared  with  unity  and  the  variance  of  H(k)  is 
small. 

°ck  “  '  C(kAx))2 

=  ~~  1°“  +  C2(kAr)3  (2.31) 

R(k)  J 

W 


The  value  of  C(0)  is  identical  with  the  previously  given 
mean- square  intensity  estimate  if  the  same  sample  mean  is 
subtracted,  and  the  previous,  result  with  low  mean  sample 
rate  and  large  N  reduces  to  the  usual  formula: 

a2  =  1 of 

c°  N  (2.32) 

2.4  POWER  SPECTRUM 

2.4.1  Computation  From  the  Autocovariance 

A 

Once  the  discrete  autocovariance  estimate  C(k) 
has  been  computed,  the  two-sided  power  spectrum  estimate 
follows  using  a  discrete  Fourier  transform  [4,7)  as 

S(iAf)  *  At  [ C (0)  +  2J1  w(k)  C(k)  costijfc1)) 

ksl  M 


(2.33) 
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where 


M  = 


Tmax 

At 


Af 


1 

2MAt 


1 

2t 

max 


f  =  iAf 


and  w(k)  is  a  smoothing  window  function  which  is  unity 
at  k=0  and  which  decreases  to  zero  at  k=M.  The  discrete 
transform  may  be  implemented  to  a  fast  Fourier  transform 
with  certain  precautions  regarding  the  addition  of  zeros 
as  discussed  by  Brumbach  [7] . 

Several  observations  are  appropriate.  Details 
concerning  these  observations  may  be  found  in  the  literature 
[2,4,7,9,17]  since  they  are  common  to  conventional  digital 

A 

spectrum  estimation.  First,  S(f)  is  an  even  periodic  func¬ 
tion  with  period  1/At.  Since  the  estimate  replicates 
itself  for  f>lAr,  the  usual  Nyquist  criterion  is  a  neces¬ 
sary  condition  on  At  to  avoid  aliasing  distortion,  just 
as  it  is  in  periodic  sampling.  The  additional  constraint 
on  At  mentioned  in  section  2.3.2  is  discussed  further 
below. 

Optimum  selection  of  the  shape  of  the  window 
function  is  an  art  involving  knowledge  of  the  spectrum 
shape  (either  a  priori  or  by  trial  and  error  procedure) 
and  the  measurement  objectives.  For  smooth  broadband 
turbulence  spectra  without  spikes  due  to  periodic  com¬ 
ponents  the  simple  Bartlett  window  given  by 

w(k)  -  1  -  -Jj-  ,  k  <  M  (2.34) 

=  0  otherwise 
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should  be  adequate.  The  effect  of  the  window  is  a  running 
weighted  average,  or  convolution,  in  the  frequency  domain 
which  reduces  the  error  variation  of  the  estimate.  At 
the  same  time,  bias  error  is  introduced  in  the  form  of  loss 
of  spectral  resolution  and  the  expected  value  of  the 

A 

estimate,  assuming  C(k)  is  unbiased,  is 


S(f)  =  S(f)*W(f) 


(2.35) 


where  *  denotes  convolution  and  W(f)  is  the  Fourier  trans¬ 
form  of  w(t) .  For  the  Bartlett  window  given  in  (2.34) 


K(f)  .  MAT  [  ]  2 


(2.36) 


the  "frequency  resolution",  R,  may  be  taken  as  1/MAt,  the 
spread  between  the  peak  and  the  first  zero  of  W(f) . 

Several  other  window  functions  having  better 
resolution  and/or  side- lobe  characteristics  are  described 
in  the  literature.  The  "Hamming"  window  is  a  common  com¬ 
promise  favorite  which  performs  better  than  the  Bartlett 
window  in  most  cases.  It  is  given  by 

W(k)  =  0.54  +  0.46  cos  (j-p) ,  k  <  M  (2 

=0  otherwise 
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A  final  observation  concerns  the  temptation  on 
the  part  of  some  to  compute  S(f)  at  frequency  intervals 
less  than  1/(2At)  because  the  result  is  smoother  and  more 
pleasing  looking  than  the  somewhat  jagged  appearance  of  the 
result  of  using  (2.33).  Giving  into  this  temptation  does 
not  produce  more  useful  information  or  resolution  since  the 
width  of  W(f)  generally  exceeds  l/(2MAx).  On  the  other  hand, 
natural  spectra  are  usually  ‘'smooth",  and  the  jagged  appear¬ 
ance  resulting  from  errors  is  not  disguised  when  the  frequency 
spacing  is  maintained  large  enough  for  the  estimates  to  re¬ 
main  nearly  independent. 


2.4.2  Bias  Effect  of  At  Selection 


Under  the  conditions  AAt  <<  1  and  N  >>  MXAx, 
result  (2.26)  becomes 


c  (k)  -  £  C(T)  ji  Reel  (54JII-)  dx 


(2.38) 


where 


Rect  (t)  =  1,  -i/2  <  t  <  1/2 
*  0  otherwise 


(2.39) 


Equation  (2.38)  is  a  convolution  integral  which  produces 
the  effect  in  the  spectrum  estimate  of  multiplication  by 
the  transform  of  the  Rect  function.  Manipulation  with 
Fourier  transform  relations  yields 


S(f)  -  l 

-  CO 


S(f-  ^)  Sin  [TrATtf-g7)] 


(2.40) 
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Equation  (2.40)  reveals  the  usual  nature  of  aliasing  error 
wherein  components  of  S(f)>l/2AT  fold  and  add  to  the  esti¬ 
mate.  The  equation  also  shows  the  attenuation  of  all 
spectral  components  by  a  (sin  x)/x  function  which  is  down 
to  a  factor  of  0.637  at  the  1/2At  folding  frequency. 

2.4.3  Effect  of  Non- Zero  Mean  and  Linear  Trends 


The  finite  resolution  imposed  by  the  finite 
width  of  the  spectral  window  function  (M<°°)  means  that  any 
discrete  frequency  component  in  the  data,  which  ideally 
would  produce  an  impulse  in  the  spectrum,  will  result  in 
a  replication  of  the  spectral  window  at  the  location  of 
the  frequency  spike  with  bias  distortion  of  nearby  esti¬ 
mates.  Even  in  the  absence  of  periodic  velocity  fluctua¬ 
tions,  the  steady  mean  velocity  component  has  the  same 
effect  at  the  frequency  origin.  The  presence  of  a  large 
linear  trend  in  the  data  which  results  from  a  sinusoidal 
component  too  low  in  frequency  to  be  resolved  also  has 
the  same  effect.  For  this  reason,  it  is  common  spectral 
estimation  practice  to  remove  both  the  sample  mean  from 
the  data  and  to  further  correct  for  linear  trends  [4] . 

The  necessity  for  linear  trend  removal  in  tur¬ 
bulence  measurements  does  not  seem  apparent,  and  the  theory 
concerning  such  a  step  in  the  case  of  random  sampling  has 
not  been  considered.  In  addition,  the  exact  nature  of  the 
effect  of  removing  the  sample  mean  is  also  slightly 
unclear  when  short  data  segments  are  used.  We  believe 
that,  when  the  data  segments  are  sufficiently  long  to 
satisfy  other  requirements,  any  low  frequency  bias  errors 
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resulting  from  subtraction  of  the  Sample  mean  instead 
of  the  true  mean  will  be  confined  primarily  to  the  zero 
frequency  estimate  with  some  leakage  to  the  first  fre¬ 
quency  estimate  (f=Af)  through  the  greater  than  l/(2MAx) 
width  of  the  spectral  window. 


2.4.4  Variance  of  the  Spectral  Estimate 


In  Appendix  II  we  derive  the  variance  of  the 
spectrum  estimate  under  the  assumption  that  X/2B  <<  XMAx  <  1. 
The  derivation  also  assumes  that  At  is  small  enough,  M  and 
N  are  large  enough,  and  the  true  mean  has  been  subtracted 
so  that  the  resulting  estimate  is  unbiased  and  the  histo¬ 
gram  H(k)  is  approximately  a  constant  given  by  NXAx .  The 
result  is 


a 


2 

f 


(S(iAf)-S(iAf))2 

M 

-  jjj1  I'  «*Ck)  [a1*  +  C^kAt)] 


(2.41) 
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This  result  shows  the  mean  square  error  has  a  maxima  at 
the  zero  frequency  and  the  folding  frequency  with  values 
which  decrease  to  a  value  approximately  one  half  as  large 
for  frequencies  away  from  i=(0,M). 

If  the  summation  is  viewed  as  a  numerical  inte¬ 
gration,  the  variance  away  from  the  frequency  origin  may 
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be  evaluated  approximately  by  neglecting  the  small 
contribution  from  the  C2(kAx)  term  as 

,4 


a,  - 


„  2a 


f 


w 


(T)dx 


(2.42) 


which  for  the  Bartlett  window  gives 

-2  „  2o4MAt 
af  "  3NX - 


(2.43) 


This  result  shows  clearly  the  price  paid  for  increasing 
MAt  to  obtain  better  spectral  resolution. 

We  may  compare  our  error  prediction  result  with 
that  given  by  Jones  [17]  for  spectrum  estimation  with 
missing  observations.  Putting  his  results  in  the  same 
notation  and  using  our  result  for  the  value  of  the  histo¬ 
gram  H(k)  we  obtain 


f j  NAAt 


M 

I 

o 


w2  (k) 


(2.44) 


where  replacement  by  a  continuous  integral  gives 


2S2(f) 

NAAx2 


MAt 


w2 (x)dx 


(2.45) 


This  result  agrees  exactly  with  ours  under  the  condition 
that  the  true  spectrum  is  white  so  that  o2  =  S(f)/Ax.  Our 
result,  however,  indicates  that  the  error  level  is  not 
proportional  as  a  function  of  frequency  to  the  true  spec¬ 
trum  as  Jones’  result  implies. 

The  usual  practice  of  indicating  a  constant 
confidence  interval  on  a  logarithmic  spectrum  plot  is 
obviously  not  appropriate  when  the  error  level  is  not 
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proportional  to  the  spectrum  itself.  It  now  becomes 
more  appropriate  to  indicate  a  constant  confidence  inter¬ 
val  on  a  linear  spectrum  plot.  Supporting  evidence  for 
this  statement  in  addition  to  our  own  may  be  found  in 
careful  observation  of  the  digital  simulation  results  in 
Jones'  previously  referenced  paper.  Unfortunately  our 
analysis  is  supported  by  experiment  with  rather  grim 
consequences  for  the  estimation  of  low  level  spectral 
tails. 


2 v 4 . 5  Post  Estimate  Smoothing 

Nature  has  provided  some  assistance  in  the 
dilemma  which  results  when  the  constraints  of  resolution, 
error,  and  measurement  time  are  all  considered.  For 
many  phenomena,  a  constant  resolution  on  a  log- frequency 
plot  is  more  meaningful  than  is  the  constant  resolution 
which  results  in  the  preceding  discussion.  Various  compu¬ 
tation-saving  methods  of  achieving  resolution  at  least 
appoximately  proportional  to  frequency  which  reduce  error 
variation  at  high  frequencies  have  been  devised  for  con¬ 
ventional  spectral  analysis  [4,24]. 

Computationally  efficient  methods  of  achieving 
non-uniform  spatial  resolution,  appropriate  for  random 
sampling  should  be  pursued,  perhaps  by  extension  of  low- 
pass  digital  filtering  to  randomly  sampled  data  [25],  but 
in  their  absence  the  error  benefits  may  still  be  obtained 
by  additional  post- estimation  smoothing.  This  is  accom¬ 
plished  by  replacing  each  spectral  estimate  with  an 
average  over  a  number  of  the  original  neighboring  estimates 


25 


AEDC-T  Ft -74-53 


proportional  to  the  frequency  of  the  estimate.  For 
example,  if  Af  is  100  Hz,  the  resolution  is  200  Hz,  and 
there  are  500  estimates  from  100  Hz  to  50  KHz,  it  might 
be  appropriate  to  retain  no  better  than  10%  resolution. 

This  would  result  in  leaving  the  first  10  estimates 
unchanged,  averaging  ten  estimates  on  either  side  along 
with  the  10  KHz  estimate,  and  averaging  45  estimates  on 
either  side  with  the  45  KHz  estimate. 

If  all  of  the  spectral  estimates  were  uncorrela¬ 
ted,  the  variance  of  the  post- smoothed  estimate  would  be 
reduced  by  the  inverse  of  the  number  of  estimates  averaged. 
Since  the  resolution  generally  exceeds  the  Af  spacing, 
the  actual  reduction  in  variance  is  less  due  to  correla¬ 
tion  between  adjacent  estimates.  This  effect  has  not 
been  analyzed  in  detail,  but  a  factor  of  2  seems  conser¬ 
vative  when  the  shape  of  the  Bartlett  window  is  considered. 

2.5  PROCESSOR  DEAD  TIME  EFFECT 

Our  idealized  assumptions  do  not  provide  for 

the  significant  fact  that  all  LV  burst- counter  processors 

exhibit  dead  time,  d,  during  which  the  processors  make 

the  "instantaneous"  sample  measurement  and  verifies  the 

validity  of  the  measurement  with  various  data  checks  which 

have  been  devised.  The  significance  lies  in  the  fact  that 
2  A 

the  minimum  value  of  At  is  2d  if  CC1)  is  not  to  be  biased. 

It  would  thus  seem  that  the  electronic  signal  processor 
dead  time  limits  the  highest  spectral  frequencies  which 
may  be  observed. 

2  This  constraint  could  be  reduced  slightly  by  restric¬ 
ting  the  band  of  delays  to  less  than  ±At/2  about  each  value 
kAt  at  the  expense  of  loss  of  data  and  more  complex  data 
processing  requirements. 
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The  minimum  possible  processor  dead  time  is  bounded 
by  the  LV  system  parameters:  it  is  the  time  required  for 
passage  of  a  scattering  particle  through  the  probe  volume. 

2.6  SEGMENTED  DATA  -  NONSTATIONARITY 

The  ergodic  assumption  includes  statistical  sta- 
tionarity  of  both  the  turbulence  and  the  random  sampling 
process.  In  many  practical  situations  one  or  both  sta- 
tionarity  assumptions  will  be  violated.  The  analysis  of 
nonstationary  random  processes  is  difficult,  and  generally 
requires  characterization  models  and  parameters  which  are 
not  available.  Evans  and  McCarty  [26]  have  discussed  such 
problems  at  length  for  the  case  of  periodically  sampled 
data. 

The  problem  we  face  is  that  in  implementation 
we  must  choose  between  obtaining,  storing  and  processing 
one  large  amount  of  data  N  as  one  continuous  record;  or 
computing  many  estimates  from  short  data  segments  and 
averaging  the  results;  or  doing  something  that  is  a  compro¬ 
mise.  We  may  further  subdivide  the  choices  by  pointing 
out  that  the  results  will  be  different  if  data  segments 
are  chosen  to  be  of  equal  length  in  time,  N'At,  or  equal 
length  in  number  of  data  words  N.  Under  the  stationarity 
assumptions  our  theory  is  directly  applicable  to  individual 
segment  estimates  with  identical  numbers  of  samples  N. 

The  succeeding  averaging  of  Ns  segments  will  not  affect 
the  predicted  bias  errors,  but  it  will  reduce  variance  in 
a  manner  simply  obtained  by  replacing  N  by  NNs  in  the 
equations.  It  is  expected  that  differences  in  implementations 
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will  be  minor  when  the  random  process  are  stationary, 
and  that  approximately  correct  predictions  will  be  obtained 
by  replacing  N=AT  when  it  is  desirable  to  do  so  for  design 
purposes . 

In  the  absence  of  theoretical  guidelines,  it 
appears  advantageous  to  make  use  of  the  fact  that  in  many 
nonstationary  situations  the  statistics  may  be  ’’slowly 
varying”  in  such  a  manner  that  the  processes  are  quasi¬ 
stationary  over  short  time  segments.  Under  these  ill 
defined  conditions  the  use  of  statistical  averages  may 
provide  valid  answers  in  terms  of  the  time  dependent 
statistics,  although  the  variance  of  the  individual  esti¬ 
mates  may  be  unusably  large  and  the  effect  of  averaging 
the  individual  estimates  to  obtain  a  long-time  estimate 
is  not  known. 
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SECTION  III 

EXPERIMENTAL 


3.1  INTRODUCTION 

An  experimental  investigation  was  conducted  in 
parallel  with  theoretical  studies  to  guide  the  direction 
of  the  theoretical  research  and  to  verify  results  as  they 
were  derived.  A  hybrid  electronic  simulation  facility  was 
designed  and  constructed  to  produce  randomly  sampled  exper¬ 
imental  data  for  the  investigation  and  computer  programs 
were  written  to  process  the  data  in  records  (segments)  of 
250  samples  each. 

Several  stages  of  experimental  work  followed. 

The  first  stage  involved  a  comparison  of  the  performance 
of  different  spectrum  estimators  (see  Appendix  I)  and  the 
latter  stages  were  studies  of  a  particular  estimator  (see 
section  2.4)  with  emphasis  on  how  its  error  behaved  as  a 
function  of  average  sample  rate  and  bandwidth  as  well  as 
the  estimator  parameters. 

A  detailed  description  of  the  simulation  facility, 
the  data  processing  procedures  and  a  copy  of  the  subroutine 
RASPEC  implementing  the  spectrum  estimator  of  Section  2.4 
are  all  included  in  Appendix  III.  A  description  of  the 
main  experimental  results  follows. 
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3 . 2  RESULTS 

3.2.1  Preliminary  Results 

Initially,  three  spectrum  estimators  were  pro¬ 
grammed  and  tested  with  the  same  randomly  sampled  data. 
Brief  descriptions  of  the  first  two  estimators  and  the 
test  results  are  given  in  Appendix  I.  The  third  esti¬ 
mator,  which  gave  the  least  qualitative  error  with  less 
computational  effort,  is  described  in  Section  2.4. 

Other  experimental  work  followed  using  the 
third  estimator.  A  low  frequency  sinusoid  (200  Hz)  was 
randomly  sampled  at  two  different  average  sample  rates 
and  this  data  produced  clearly  discernable  spikes  at  the 
fundamental  frequency  with  only  250  data  samples.  Both 
bandpass  and  lowpass  random  processes,  generated  by  pass¬ 
ing  white  noise  through  controllable  known  filters,  were 
then  randomly  sampled  for  a  number  of  average  sample  rates, 
ranging  from  below  the  Nyquist  rate  to  above  the  Nyquist 
rate . 


3 . 2 . 1 . 1  Experimental  Error  as  a  Function  of  the  Number  of 
Data  Words 


The  decrease  in  experimental  spectrum  estimate 

A 

error,  e^s,  with  an  increase  in  the  number  of  data  words 
processed  is  shown  in  Figure  1,  a.  and  b.,  where 
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The  subscript  s  indicates  that  this  is  the  error  in  the 
stopband  which,  in  this  case,  is  the  band  from  500  Hz  to 
2000  Hz  where  the  true  spectrum  is  assumed  to  be  zero. 
There  are  151  estimates  in  this  band  at  10  Hz  intervals. 

A 

Qualitatively,  elg  is  the  square  root  of  the  average  error 
squared,  normalized  to  the  peak  of  the  spectrum  estimate. 

The  same  pair  of  estimates,  plotted  on  a  dif¬ 
ferent  scale  (figure  1,  c.  and  d.),  indicates  that  the 
estimator  produced  negative  and  positive  errors  that  are 
symmetrical  about  the  true  spectrum.  This  is  clearly  seen 
in  the  stopband. 


3.2. 1.2  Experimental  Error  as  a  Function  of  Mean  Sample  Rate 


In  the  series  of  estimates  shown  in  Figure  2,  the 
number  of  data  words  (N=248)  and  the  number  of  data  segments 
(Ns=20)  are  held  constant  and  the  mean  sample  rate,  X,  is 

A  A 

increased.  The  reduction  in  els  is  indicated,  where  els 
is  defined  just  as  in  the  previous  section.  Note  that  in 
the  main  lobe  of  the  bandpass  spectrum,  the  estimate  error 
qualitatively  appears  to  increase  for  the  higher  values  of  X 
where  a  =  X/2B  >>  1. 


3. 2. 1.3  Aliasing  Error  and  the  Effect  of  Smoothing 

Aliasing  error  is  shown  in  Figure  3.  In  Figure 
3a. ,  the  parameter  Ax  was  purposefully  selected  large  enough 
to  cause  the  spectrum  estimate  to  replicate  itself  within 
the  scale  of  the  plot.  A  further  increase  in  Ax  results 
in  the  increased  aliasing  of  Figure  3b.  The  error  reduction 
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1. 


Figure  2.  Variation  of  Spectrum  Estimates  with  Mean  Data 
Rate. 
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Figure  3.  Spectrum  Estimates  Exhibiting  Aliasing. 
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effect  of  decreasing  the  width  of  the  window  w(k)  is 
shown  in  Appendix  I,  Figure  1-4. 

3.2.2  Parametric  Variability  Study 

One  of  the  main  theoretical  results  was  the 
derivation  of  an  expression  for  the  variance  of  the  spec¬ 
trum  estimator  of  Section  2.4.  Accordingly,  an  experi¬ 
mental  effort  was  undertaken  to  verify  this  expression 
(equation  (2.43))  by  plotting  theoretical  error  and 

experimental  error  as  a  function  of  the  mean  rate  parameter, 
X 

a  =  TB  * 

To  obtain  a  range  of  values  for  a,  random  process¬ 
es  with  two  different  lowpass  spectra  were  used,  one  with 
B=100  Hz  and  one  with  B  =  800  Hz  (nominal  values).  Four 
different  average  sampling  rates  were  used  for  each  of  the 
two  spectra  to  give  eight  values  of  the  parameter  a.  This 
experimental  data  was  then  used  to  compute  a  Mhigh  accuracy” 
spectrum  estimate  for  each  bandwidth.  For  this,  four  20 
data  segment  spectrum  estimates  were  computed  for  each  value 
of  X. -  These  estimates  were  averaged  at  each  point  in  fre¬ 
quency  to  form  an  averaged  estimate  for  each  value  of  X. 
These  four  averaged  estimates  (from  the  four  different 
values  of  X)  were  again  averaged  at  each  point  in  frequency 
to  form  the  overall  "high  accuracy"  estimate  for  that  band¬ 
width.  The  "high  accuracy"  estimates  were  then  fitted  with 
the  ideal  Butterworth  spectrum  (fourth  order)  for  each  band¬ 
width  using  a  least-mean  squares  program  to  determine  the 
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level  o£  the  Butterworth  curve  in  the  passband.  Wave 
analyzer  data  showed  that  frequencies  below  about  20  Hz 
were  attenuated  in  the  actual  spectrum  and  so  it  was  not 
a  true  Butterworth  lowpass  spectrum.  Consequently,  the 
first  six  points  of  the  estimates  were  not  used  in  the 
least-mean  squares  program.  The  plots  of  Figure  4  show 
the  "high  accuracy"  estimates  and  least-mean  square  fit 
Butterworth  curves  for  each*  bandwidth. 

The  next  step  involved  was  to  compute  experimental 
error  by  computing  error  between  a  number  of  spectral 
estimates  and  the  Butterworth  curves  described  above.  For 
each  value  of  a,  80  data  segments  were  used  in  computing 
estimates;  one  estimate  from  the  first  20  data  segments, 
and  so  forth  for  a  total  of  four  estimates.  This  resulted 
in  32  spectral  estimates.  The  experimental  error  was 
defined  over  two  frequency  bands  for  each  bandwidth.  The 
bands  were: 

Passbands 

24  Hz  <  f  <  104  Hz  100  Hz  filter 

204  Hz  <  f  <  800  Hz  800  Hz  filter 

Stopbands 

204  Hz  <  f  <  00  Hz  100  Hz  filter 

1632  Hz  <  <  6400  Hz  800  Hz  filter 
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Analytical  Butterworth  Spectra  and  "High 
Accuracy"  Spectrum  Estimates. 
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where  the  upper  and  lower  limits  were  determined  by  counting 
off  a  discrete  number  of  points  of  frequency  spacing  Af. 

The  general  expression  for  the  experimental  error 


in  the  passband. 


was 


n+m-1 


lp 


tii  ■  hi 

L  .1 

n.  i=m 


(S(iAf)  -  S(iA.f)  ) 


2] 


6(0) 


(3.2) 


where  S(iAf)  is  the  value  of  a  spectrum  estimate  at  frequency 
f  =  iAf. 

For  the  100  Hz  filter, 

/v 

S(iAf)  =  value  of  the  analytical  model  (100  Hz 
filter)  at  frequency  f  =  iAf. 

S(0)  =  value  of  analytical  model  (100  Hz  filter) 
in  the  low  frequency  portion  of  the  pass- 
band 

n  *  number  of  points  in  the  passband  over 
which  the  summation  is  taken  =  20 
m  =  first  point  in  the  passband  =  6 
Similar  definitions  apply  for  the  experimental  error  in  the 
passband  of  the  800  Hz  filter. 

The  general  expression  for  the  experimental  error 

A 

in  the  stopband,  ,  is  defined  similarly  to  the  experimental 
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error  in  the  passband,  except  the  summations  are  taken 
over  n  =  149  and  m  =  51  in  the  stopbands  of  both  filters. 

By  applying  these  definitions  to  the  32  spectrum 

A 

estimates,  two  plots  were  generated,  one  for  e,_  vs  a,  the 
normalized  mean  rate  parameter,  (Figure  5)  and  one  for  e^s 
vs  a  (Figure  6) .  For  each  value  of  a,  four  values  of  e^p 
were  computed,  one  from  each  of  the  four  spectrum  estimates. 
These  were  then  averaged  to  form  e,  and  this  is 

X  P  "*  flV  A 

plotted  as  a  solid  line  in  Figure  5.  The  values  of  e^p 
from  the  first  20  records,  and  the  second  20  records, 

e^p_2»  are  shown  as  dashed  lines  in  Figure  5  to  show  typi¬ 
cal  variability  of  the  values  about  their  average. 

To  compare  theory  and  experimental  results,  equa¬ 
tion  (2.43)  was  plotted  in  normalized  rms  form: 


e 


IP 


(3.3) 


where  a 2  =  2BS(fpea^)  has  been  substituted  using  the  nor¬ 
malized  mean  rate- parameter  a=  jg-  .  The  values  of  the 
parameters  used  in  the  experiment  were  substituted  in  the 
equation  for  computations. 

A 

The  description' of  Figure  6,  the  plot  of  elg  vs 
a,  is  essentially  the  same  as  for  Figure  5  but  the  results 
are  somewhat  different.  The  plot  shows  much  less  varia- 

A  A 

bility  between  the  empirical  average  e,  and  e-  ,  or 
a  xs “ av  l s  ~  x 

els-2  an<^  a^so  shows  much  better  agreement  between  els  (the 
theoretical  curve)  and  eis_av-  The  agreement  between  the 
empirical  and  theoretical  is  good  even  for  relatively 
large  values  of  a. 
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Figure  6.  Spectrum  Estimate  Error  in  the  Stopband 
Theoretical  and  Experimental. 


AEDC-TR -74-53 


3.2.3  Parametric  Study  of  Finite  Time  Correlation  Effects 
on  Mean  and  Turbulence  Intensity  Variability 

Experiments  similar  to  those  described  in  Section 
3.2.2  were  conducted  to  determine  whether  or  not  the  increase 
of  varability  due  to  finite  time  correlation  effects  were 
important  for  fixed  N.  The  results  are  shown  in  Figures  7 
and  8  where  experimental  measures  of  the  normalized  RMS 
deviations  of  the  mean  and  turbulence  intensity  estimates 
are  plotted  versus  a  for  N=248  and  Ns=400.  The  theoretical 
curves  are  plots  of  equation  (2.14)  (mean  estimate)  and 
equation  (2.21)  (turbulence  intensity  estimate)  in  normal¬ 
ized  form.  It  is  seen  that  estimate  error  did  increase 
for  a  >  1  but  was  less  than  predicted.  The  theoretical 
equations  appear  conservative  and  further  work  is  required 
for  exact  results. 

3.2.4  Histogram  Error 

Figure  9  shows  the  finite  time  effects  on  the 
histogram  H(k)  of  lag  products  as  a  increases.  Figure  9a 
shows  the  histogram  H(k)  and  its  expected  value,  equation 

(2.27) ,  for  comparison.  Figure  9b  shows  the  histogram 
H(k)  and  equation  (2.27)  under  the  condition  of  equation 

(2.28)  . 
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SECTION  IV 
PRELIMINARY  DESIGN 


4.1  INTRODUCTION 

In  this  section  we  provide  some  extension  and  fur¬ 
ther  interpretation  of  the  theoretical  results  of  Section  II 
in  the  form  of  a  preliminary  systems  design  guidelines  and 
provide  two  illustrative  examples.  The  discussion  is  restrict- 
ted  to  the  requirements  of  spectral  analysis  since  implemen¬ 
tation  of  mean  and  turbulence  intensity  estimates  is  straight¬ 
forward  under  the  assumptions  of  lack  of  correlation  between 
X  and  U(t) . 

4.2  PRELIMINARY  DESIGN  PROCEDURE 

There  are  many  interrelated  system  variables 
whose  proper  relationship  must  be  maintained  to  achieve 
results  which  satisfy  error  objectives  while  minimizing 
data  collection  and  processing  time.  The  same  procedures 
must  be  used  in  system  design,  experiment  planning,  and 
final  data  analysis.  The  procedure  is  applied  iteratively 
because  of  the  many  tradeoffs  and  the  general  fact  that 
a  priori  knowledge  of  the  spectrum  to  be  estimated  is 
essential  in  the  analysis.  With  the  exception  of  the 
alteration  of  the  relationships  due  to  random  sampling,  the 
design  process  is  similar  to  that  involved  in  conventional 
spectral  estimation.  Detailed  design  discussion  for  con¬ 
ventional  estimation  has  been  provided  by  Blackman  and  Tukey 
[2]  and  Parzen  [28]  . 


46 


AEDC-TR -74-53 


4.2.1  Bandwidth  -  Selection  of  At 


Naturally  produced  spectra  are  not  sharply  band- 
limited  and  some  aliasing  bias  error  will  generally  have 
to  be  tolerated.  When  designing  a  system  for  turbulence 
spectrum  estimation,  electrical  engineers  accustomed  to 
the  use  of  3  db  bandwidth  must  exercise  caution  and  communi¬ 
cate  clearly  with  the  aerodynamicists  about  the  desired 
high  frequency  response.  By  way  of  example  we  note  that 
a  simple  one-pole  lowpass  spectrum  with  3  db  bandwidth  of 
2.5  KHz  and  log- spectrum  tail  slope  of  -2  has  an  equivalent 
power  bandwidth  of  3.92  KHz  and  is  only  26  db  down  at  50  KHz. 
For  this  example,  the  factor  C2  used  in  4.2.4  below  is  0.0785. 

The  aliasing  error  equation  (2.40)  indicates  that 
aareful  consideration  must  be  given  to  both  the  usual  addi- 

j 

tive  aliasing  error  and  the  multiplicative  error  due  to 
random  sampling.  In  a  typical  case,  the  maximum  error  occurs 
at  a  frequency  less  than  the  folding  frequency  1/2At  because 
the  multiplicative  error  and  the  additive  errors  may  cancel 
at  the  folding  frequency. 

4.2.2  Resolution,  Frequency  Spacing,  and  Number  of  Estimates 

The  frequency  resolution  R  is  an  ambiguous  quantity 
related  to  window  shape,  method  of  definition,  and  the  lack 
or  presence  of  discrete  frequency  spikes.  For  broadband 
spectra  with  R  defined  as  the  half  width  of  the  spectral 
window  measured  to  the  first  zero  and  with  Af  selected  as 
1/2MAt  where  M  is  the  maximum  lag  number  and  also  the  number 
of  spectral  estimates,  then  R=2Af  when  a  Bartlett  or  Hamming 
window  is  selected. 
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The  value  of  M  is  given  by 

M  =  2Af  At  (4.1) 

with  typical  values  being  256,  512,  1024.  It  is  important 
to  note  that  it  is  a  simple  matter  to  reduce  the  value  of 
M  after  the  autocovariance  estimate  has  been  made  if  it 
proves  desirable  to  sacrifice  resolution  for  reduced  varia¬ 
bility  error.  In  other  circumstances,  post-estimate  smooth¬ 
ing  may  prove  desirable.  However,  it  M  is  chosen  inadequately 
small,  no  simple  remedy  is  available. 

4.2.3  Data  Segment  Length 

In  some  implementations  it  is  desirable  to  process 
data  in  segments  and  subtract  the  sample  mean  for  the  segment 
from  each  sample.  If  the  segment  length,  N,  is  too  small  the 
turbulence  intensity  estimate  may  be  biased  low  according 
to  equation  (2.1-6.).  The  bias  effects  on  the  spectrum  are 
unknown,  but  are  believed  to  be  confined  to  the  zero  fre¬ 
quency  estimate  with  leakage  through  the  window  function  to 
the  estimate  at  f=At. 

An  additional  constraint  on  N  is  imposed  by  equation 

i 

(2.28)  to  prevent -excessive  variance  of  the  spectral  esti¬ 
mate  due  to  low  values  of  the  histogram  H(k)  for  the  larger 
values  of  k.  This  constraint  is  generally  satisfied  by 

N  >  M  (4.2) 

for  sample  rates  X/2B  <  1,  AAt  <<  1.  Note  that  this  con¬ 
straint  may  be  violated  after  system  implementation  by  an 
over-eager  experimenter  who  increases  \  to  decrease  varia¬ 
bility  error.  The  inclusion  of  a  comparison  of  H(l)  and 
H(M-l)  as  a  check  is  recommended. 
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4.2.4  Variability  Error 


Once  bias  error  effects  are  included  properly  in 
the  system  design  and  experimental  execution,  the  primary 
concern  is  the  statistical  error  due  to  finite  data  collec¬ 
tion.  Equation  (2.42)  may  be  expressed  logarithmically  for 
worst  case  error  analysis  in  terms  of  the  normalized  RMS  spec¬ 
trum  error 


e  = 


S<W 


(4.3) 
The  worst 


where  S(fpea^)  is  the  peak  spectrum  value, 
case  (largest  value  of)  equivalent  power  bandwidth  B 


B 


is  used  to  obtain 


(4.4) 


log  Ej  =  %  log  (2C1C2)  -  %  log  a  -  fc  log  (B)  (4.5) 

where 

r  =  _L_rMAT  a  =  ^4 

'■‘1  MAt^  o  w2(t)dr  2B 

C-  =  2AtB  g  =  Ml 

2  P  M 


This  result  is  plotted  in  Figures  10  and  11  in  db  using 
typical  values  of  =  1/3  (Bartlett  window)  and  C2  =  0.0785. 

Equation  (2.42)  may  also  be  expressed  in  terms  of 
the  total  data  collection  time  T  (not  including  any  process¬ 
ing  time  between  data  segments)  by  replacing  NNs  with  the 
expected  number  of  samples  XT: 

log  Ej  =  \  log  (2C ^ )  -  log  (a)  -  \  log  (y)  (‘4'6') 
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Normalized  Number  of  Data  Samples,  $ 

Figure  10.  Normalized  RMS  Spectrum  Error  o  vs 
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Normalized  Number  of  Data  Samples,  B 
Figure  11.  Normalized  RMS  Spectrum  Error  vs  6. 
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where  and  a  are  the  same  as  in  (3.6)  and 

T 

Y  'MAt  (4.7) 

The  result  is  plotted  in  Figures  12  and  13. 

We  examine  the  significance  of  the  error  equations 
in  examples  in  Section  4.3  below.  First  there  are  two  more 
primary  considerations. 

4.2.5  Amplitude  and  Time  Error 

The  present  analysis  has  not  included  the  effects 
of  amplitude  errors  due  to  quantum  and  electronic  noise, 
velocity  gradients  and  particle  lag  effects,  and  timing 
errors  in  both  the  time-between-samples  and  the  signal 
period  counters.  In  addition,  fast  hardware  reciprocation 
(period  to  scaled  velocity)  conversion  has  been  assumed 
without  regard  to  error.  None  of  these  effects  were  observed 
in  the  experimental  results  reported  herein  because  they 
were  negligible  or  non-existent  in  the  simulation.  However, 
they  should  be  expected  in  LV  measurements  of  high  velocity 
gas  flow.  Further  discussion  is  presented  in  Section  5. 

4.2.6  Computation  Time  and  Memory  Considerations 

The  theory  section  does  not  concern  itself  with 
the  time  required  to  compute  the  autocovariance  estimate. 

The  time  required  to  compute  the  required  discrete  Fourier 
transform  to  obtain  the  spectral  estimate  is  small  and  is 
of  secondary  importance.  The  critical  parameter  is  the  rate 
^  at  which  lag  products  must  be  formed,  indexed,  and 
accumulated.  This  rate  parameter  determines  what  form  of 
implementation  is  possible  in  terms  of  avaliable  software 
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Normalized  Data  Collection  Time,  a 
Figure  12.  Normalized  RMS  Spectrum  Error  a  vs  Y . 
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100  1000  10,000  100,000  1,000,000 


Normalized  Data  Collection  Time,  a 
Figure  13.  Normalized  RMS  Spectrum  Error,  e  vs  y. 


AEOC-TR-74-S3 


AEDC -TR-74-53 


and  hardware. 

The  factors  which  determine  are  the  mean  data 

rate  A  and  the  maximum  frequency  resolution  selection  of 

MAx.  The  average  lag  product  rate  is  obtained  simply  as 

the  product  of  A  and  the  expected  number  of  samples  AMAx 

available  at  any  instant  in  the  preceding  MAx  seconds. 

2 

Adding  the  rate  of  u^  terms  gives 


1 1>  =  A2MAt  +  A  (4.8) 

For  segmented  data  implementations,  the  average  number  of 
product  computations  per  data  segment  Nip  may  be  obtained 
from  the  expected  value  of  the  histogram  by  summation: 


=  N  +  NAAt 


M 

T 


(AAt) 


(4.9) 

The  ( AAx ) 2  term  is  small  in  (3.4)  when  the  suggested 
condition  AAx  <<  1  is  satisfied.  Division  of  (3.9)  by  N/A , 
the  average  time  required  to  obtain  one  data  segment,  pro¬ 
duces  equation  (3.8)  except  for  the  (AAx)2  term  which 
results  from  the  loss  of  available  lag  products  for  which 
the  time  pairs  cross  a  segment  boundary. 

The  system  memory  requirements  vary  considerably 
with  implementation.  The  minimum  possible  data  memory  occurs 
when  only  the  last  JMAX  samples  and  occurrance  times 

(IK,  Ui-1’  ***  ^i-JMAX^  are  reta^ne<^  *n  real  time  processing. 
The  quantity  JMAX  is  also  useful  in  segmented  data  implemen¬ 
tations.  This  number  is  chosen  statistically  to  insure  that 
most  of  the  sample  pairs  with  time  differences  less  than 
Mat  will  be  retained  while  most  of  those  with  excessive 
time  differences  will  not  even  be  tested.  This  benefit 
is  obtained  by  processing  lag  product  pairs  along  diagonals 
of  the  segment  product  matrix  (IK  Ui+j ,  IK+1  Ui+i+j »  etc.) 
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and  stopping  when  JMAX  is  obtained  in  the  second  subscript. 

The  value  of  JMAX  is  obtained  by  noting  that  for 
a  Poisson  sampling  process  the  expected  number  of  samples 
which  will  occur  in  time  MAt  is  XMAt,  and  the  variance  of 
this  number  is  also  XMAt,  for  values  of  JMAX  >>  1,  the 
Poisson  probability  law  has  an  envelope  which  is  nearly 
Gaussian  and  the  ratio  of  the  RMS  deviation  to  the  mean 
becomes  small.  The  mathematics  are  far  from  exact,  but  we 
may  use  the  usual  3a  law  for  Gaussian  statistics  to  obtain 
the  following  condition  which  assures  that  very  few  valid 
product  pairs  are  missed: 

k 

JMAX  =  XMAt  +  3  (XMAt)^  (4.10) 

When  XMAt  >  9,  equation  (4.10)  gives  a  value  of  JMAX  less 
than  2xMAt. 

4.3  EXAMPLES 

4.3.1  Problem  Specification 

The  following  example  is  artificial  but  illustra¬ 
tion.  We  assume  that  the  spectrum  to  be  measured  will  be 
bounded  from  above  in  the  high  frequency  portion  by 

Sl(£)  =  T  wf  ,2  ’  {  >  *0  ’  2'5  KHz  '  C4.ll) 

■  o 

The  LV  signal  process  is  assumed  to  provide  accurate 
measurements  of  the  velocity  samples  and  time  differences  with 
a  system  dead  time  of  2.5  ysec  maximum.  The  measurement  objec¬ 
tives  are  given  as: 

a)  determine  details  of  S(f)  to  50  KHz 
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b)  resolution  =  100  Hz 

c)  data  collection  time  =  30  seconds 

d)  use  on-line  real-time  minicomputer  program 
without  external  high  speed  memory 

6)  provide  maximum  error  at  each  frequency  less 
than  5%  of  the  true  spectrum 

f)  use  mean  data  rate  less  than  1000  samples/sec. 

It  is  impossible  for  any  system  to  be  designed  which 
will  meet  the  measurement  objective  specifications  since  the 
estimate  error  is  statistical  and  cannot  be  specified  by  a 
maximum  value.  The  usual  spectrum  estimation  practice  of 
specifying  confidence  interval  using  chi-squared  statistics 
is  not  appropriate  since  the  random  estimate  error  is  not 
a  positive  random  variable.  At  the  present  time  only  the 
RMS  error  may  be  specified.  With  the  change  of  e)  to  "pro¬ 
vide  bias  error  of  less  than  51  and  RMS  error  of  2%  of 
the  true  spectrum  at  each  point"  the  specifications  become 
legitimate  in  form  but  unattainable.  This  may  be  verified 
by  applying  the  previously  provided  equations  and  guidelines. 

The  next  step  is  to  compromise  on  the  specification. 
In  what  follows  we  arbitrarily  relax  some  of  the  constraints 
to  obtain  two  more  realistic  examples.  The  results  are  arbi¬ 
trary  design  selections  which  cannot  be  further  optimized 
without  priority  design  objective  tradeoff  specifications i 
The  difference  in  the  results  for  the  two  examples  clearly 
illustrates  the  need  for  priority  assignment. 

4.3.2  Low  Frequency,  High  Resolution 

The  specifications  are  to  be  applied  only  to  the 
frequency  range  100  Hz  to  2500  Hz.  The  measurement  time 

f 
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objective  of  30  seconds  is  to  be  relaxed  if  necessary.  It 
is  assumed  that  no  strong  single -frequency  spikes  or  steep 
steps  are  present  which  would  provide  excessive  bias  error 
at  near-by  estimates  through  the  convolution  of  the  spectral 
window  function  with  the  true  spectrum. 

The  first  selection  of  parameters  is  given  in 
Table  1  under  "Trial  One".  In  the  first  trial  all  objectives 
are  met  except  for  T.  The  selection  of  e-^  as  23  db  is  17  db 
below  the  2500  Hz  3  db  value  with  an  additional  1.5  db  assumed 
for  the  unknown  shape  of  the  low-frequency  portion  of  the 
spectrum  and  1.5  db  for  experimentally  indicated  difference 
of  actual  error  and  theoretical  error  in  the  passband.  The 
folding  frequency  MAf  is  sufficiently  high  that  the  additive 
aliasing  error  from  the  17.5  KHz  portion  of  the  assumed 
tail  is  4%..  This  is  adequate  even  if  the  predicted  2.5%  nega¬ 
tive  error  e~  from  random  sampling  effects  does  not  materialize. 
Unfortunately  the  required  y  and  T  are  too  large  for  consideration. 

Table  1  also  presents  a  second  iteration  in  which 
the  resolution  requirement  is  reduced  a  factor  of  two  by  de¬ 
creasing  At  and  increasing  the  folding  frequency.  This  change 
makes  z^  negligible  but  does  not  reduce  T  by  a  sufficient 
amount.  The  second  trial  also  shows  the  effect  of  increasing 
the  data  rate  to  7.85  K/sec  to  give  a  =  1 .  Experiment  indicates 
that  there  would  be  little  or  no  advantage  in  increasing  X 
beyond  the  a  =  1  value  for  passband  error.  Unfortunately, 
the  mean  lag  product  rate  ^  is  too  high  for  continuous  online 
processing  but  other  alternatives  requiring  only  a  digital 
tape  recorder  are  possible  as  described  in  Section  V. 


58 


AEDC -TR-74-53 


Example  1 


Example  2 


Trial  One 

Trial  Two 

Af (Hz) 

SO 

100 

2000 

R(Hz) 

100 

200 

4000 

M 

200 

200 

50 

MAf  (Hz) 

10K 

20K 

100K 

At (sec) 

50y 

25y 

5  y 

A (/sec) 

IK 

7.85K 

23. 5K 

fQ(Hz-3db) 

2.5K 

2.5K 

2.5K 

B  (Hz) 

3.93K 

3.93K 

3-.93K 

a 

0.127 

1.0 

3.00 

'K/sec) 

11K 

316K 

138K 

1/iHsec) 

90y 

3 . 2  y 

7 . 24y. 

■i(db) 

-23 

-23 

-35.5 

•S 

4% 

negligible 

11% 

e2 

2.5% 

negligible 

-10% 

e2 

1.5% 

negligible 

11 

Y 

excessive 

26K 

900K 

T 

excessive 

130 

225 

NN 

s 

excessive 

1020K 

5288K 

Table  1.  Trial  Parameters  for  Example  Problems 
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4.3.3  High  Frequency  Spectrum  Tails 

In  order  to  obtain  2%  error  at  50  KHz  the  value 
of  relative  to  the  constant  lowpass  level  of  the  assumed 
spectrum  S1(f)  would  be  -43  db  since  S^f)  is  -26  db  at 
50  KHz.  A  few  trial  parameter  selections  shows  this  error 
level  is  very  difficult  to  reach.  As  a  demonstration. 

Example  2  provides  a  parameter  set  which  gives  -35  db  for  e j 
with  no  theory  compensation  deemed  necessary.  The  value  of 
At  =  Sysec  is  the  smallest  allowed  by  the  signal  processor 
dead  time  without  additional  bias  errors.  We  note  that  in 
this  example  the  lag  product  rate  $  is  less  than  in  the 
previous  example  because  of  the  considerable  loss  of  resolu¬ 
tion  in  the  reduction  of  MAt.  The  total  number  of  computa¬ 
tions  required  is  larger,  however.  The  bias  errors  and 
z~2  are  the  values  at  50  KHz  and  appear  to  nearly  cancel.  The 
theory  concerning  has  not  been  verified  experimentally 
and  the  value  of  e2  has  not  been  checked  at  lower  values  of 
f.  The  bias  error  could  therefore  be  as  high  as  11%  at  50 
KHz  which  would  be  of  the  same  order  of  magnitude  as  (35  db) . 
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SECTION  V 

DISCUSSION 

5.1  IMPLEMENTATION  ALTERNATIVES 

The  results  show  that  there  are  many  implementa¬ 
tion  alternatives,  depending  on  the  measurement  objectives 
and  the  values  of  X,  the  mean  data  rate,  and  ifi,  the  mean 
lag  product  rate. 

5.1.1  Batch  Processing  -  Segmented  Data 

In  many  cases  where  only  the  principal  turbulence 
power  distribution  is  of  interest  and  measurement  time  is  not 
critical,  the  simplest  implementation  will  be  an  arrangement 
similar  to  that  which  we  used  in  which  a  minicomputer  acts 
as  a  time  buffer  and  pre-processor  between  the  LV  electronics 
and  a  digital  tape  recorder.  Additional  features  could  be 
incorporated  such  as  segment  mean  computation  and  subtraction, 
segment  turbulence  intensity  estimate,  and  elimination  of 
data  points  lying  outside  3  or  4  a  deviation  (assumed  to  be 
LV  errors) .  This  mode  of  operation  would  allow  immediate 
access  to  first  order  statistics.  It  would  also  even  allow 
for  software  period- to-velocity  inversion  and  scaling  at  the 
data  rates  of  5  to  10  K. 

For  higher  values  of  X,  such  as  may  be  required 
for  tail  investigations  or  for  multiple  channel  processing, 
several  batch  mode  alternatives  are  possible.  The  fastest 
methods  would  include  either  a  large  core  or  a  high  speed 
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disc.3  The  inexpensive  method  would  be  to  simply  retain 
the  tape  drive,  use  slightly  modified  pre-processing  and 
transfer  programs,  but  let  the  computer  wait  between  seg¬ 
ments  the  additional  time  required  to  write  a  segment  on 
tape.  Possible  complications  due  to  any  low  frequency 
sampling  effect  are  not  expected  to  be  serious  due  to  the 
random  length  of  time  required  for  acquiring  a  segment. 

With  slightly  more  waiting  time  these  same  arguments  apply 
to  replacement  of  the  tape  unit  with  a  telephone  connec¬ 
tion  to  a  large  time  shared  machine. 

The  batch  processing  mode  is  recommended  at  the 
present  time  not  only  for  reasons  of  cost  and  simplicity 
but  also  because  enough  unknowns  still  remain  to  make  saving 
the  data  for  future  reevaluation  worthwhile  if  the  test  is 
an  important  one.  For  example,  after  seeing  the  results  the 
investigator'  may  wish  to  change  At  a  week  later  and  reprocess 
the  data. 

5.1.2  On  Line  Processing  -  Segmented  Data 

A 

This  alternative  is  similar  to  the  Batch  approach. 
The  only  difference  is  that  after  a  segment  is  read,  the 
processing  continues  until  the  lag  products  and  additions 
to  the  histogram  have  been  completed  before  another  segment 
is  read.  At  low  values  of  X  and  ij;,  this  approach  could  be 
made  nearly  real-time  if  two  memory  segments  are  used  with 
priority  interrupt  for  reading  and  writing  in  between  com¬ 
putations.  The  differences  between  this  approach  and  Batch 
processing  would  be  most  apparent  at  higher  data  rates  where 
i|>  is  increased  by  a  X2  factor,  particularly  in  cases  where 
multiple  channels  of  data  would  be  processed.  In  any  event, 

3 

This  is  currently  being  done  at  the  Lockheed  Georgia 
Company  - 
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the  memory  requirements  would  be  greater  because  of  the 
additional  programming.  Since  the  output  data  rate 
would  be  quite  low  in  this  approach,  the  division  of  the 
SUM  matrix  by  the  histogram,  the  multiplication  by  the 
window  function,  and  the  transform  operations  could  all 
be  done  remotely  on  a  time- shared  computer. 

5.1.3  On  Line  Processing  -  Continuous 

This  approach  is  feasible  for  sufficiently  low 
lag  product  rates.  It  could  even  be  used  in  Batch  processing 
if  no  data  is  skipped  in  the  recording  process  because  it  has 
the  advantage  of  requiring  less  data  memory  in  core.  The 
data  memory  is  indexed  in  a  circular  manner  with  the  most 
recent  sample  replacing  the  oldest.  The  length  of  the 
memory  array  need  be  no  longer  than  JMAX  as  described  in 
Section  4.2.6.  Unfortunately,  in  most  cases  the  achievable 
rate  would  be  inadequate  unless  long  times  T  are  used. 

Another  difference  in  this  approach  is  that  a 
preestimate  of  the  mean  would  be  required  for  subtraction 
from  each  sample  as  it  occurs.  This  has  the  advantage  of 
not  removing  low  frequency  components  in  the  first  frequency 
estimate  since  the  same  constant  is  subtracted  from  all  data. 
It  has  the  disadvantage  of  leaving  a  small  spike  of  unknown 
amplitude  at  zero  frequency  due  to  error  in  the  mean 
estimate  used  in  the  subtraction. 

5.1.4  Hardware 


This  is  a  class  of  alternatives  based  on  the 
second  formulation  of  the  autocovariance  estimate  given  by 
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equation  (2.25).  In  this  class  o£  alternatives  a  buffer 
containing  the  latest  sample  value  would  be  periodically 
read  with  reset  to  zero.  The  data  would  then  consist  of 
many  zeros  with  sample  values  sparsely  distributed.  The 

result  would  be  very  inefficient  with  memory  utilization. 

4 

The  alternate  formulation  also  has  differences  in  the 
required  error  analysis  as  previously  noted.  The  approach 
is  designated  "hardware”  since  it  is  conceivable  that  two 
commercial  correlators  could  be  modified  to  produce -the 
autocorrelation  of  the  sampled  data  sequence  and  the  indi¬ 
cator  sequence  separately  for  division  and  Fourier  process¬ 
ing  on  a  time  shared  computer.  The  general  approach  could 
also.be  implemented  as  software,  and  has  the  advantage  that 
time  differences  need  not  be  computed.  The  probable  dis¬ 
advantages  include  poorer  error  performance,  the  need  for 
hardware  period- to-velocity  conversion,  and  the  lack  of 
analysis  at  the  present  time  of  other  factors  which  make 
recording  the  data  worthwhile. 

5.2  RELATION  BETWEEN  THEORY  AND  EXPERIMENT 

The  experimental  results  of  Section  III  show 
excellent  agreement  with  the  theory  of  Section  II  with 
regard  to  spectral  estimate  varability  error  prediction 


4 

This  formulation  is  attributed  to  an  interpretation 
by  P.  Scott,  General  Electric  Co.,  of  the  results  from 
this  contract  presented  by  Mayo,  Shay,  and  Riter  at  the 
Oklahoma  State  University  Workshop  on  "Theory  and  Applica¬ 
tion  of  the  Laser  Doppler  Anemometer",  June  11,  1973. 
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in  the  tail  region.  That  this  agreement  persists  at  values 
of  X/2B  >  1  is  somewhat  surprising  in  view  of  the  violation 
of  the  derivation  assumptions  of  low  data  rate. 

The  experimental  results  also  show  the  theory  is 
inadequate  for  error  prediction  in  the  spectrum  passband  at 
rates  X/2B  >  1.  In  this  case,  the  experiments  substantiate 
intuitive  reasoning  based  on  information  theory  relevant  to 
periodic  sampling:  sampling  at  rates  in  excess  of  X/2B  is 
redundant;  the  only  way  to  increase  the  information  transfer 
concerning  the  random  process  being  sampled  is  to  increase  T. 
Continuation  of  this  intuitive  logic  shows  that  for  a  fixed 
number  of  samples,  increasing  X  above  2B  is  detrimental 
becauses  it  reduces  T.  The  quantitative  experimental  results 
of  Figure  5  are  inadquate  to  prove  that  ^  actually  exhibits 
a  minimum  at  X/2B  =  1,  but  the  possibility  is  certainly  not 
eliminated.  This  possibility  is  qualitatively  substanti¬ 
ated  by  the  sequence  of  Figure  2  where  the  values  of  X/2B 
attained  were  greater. 

We  conclude  from  our  results  that  spectrum  esti¬ 
mation  with  mean  sample  rates  much  less  than  twice  the  highest 
frequency  are  not  only  possible  but  also  desirable  when  we 
are  concerned  with  the  primary  turbulence  power  distribution. 
(The  quantity  B  is  generally  much  less  than  the  highest  sig¬ 
nal  frequency) .  On  the  other  hand,  attainment  of  usably  low 
estimate  error  in  the  spectrum  tails  will  generally  require 
many  compromises  in  terms  of  measurement  time,  resolution, 
the  use  of  attainable  but  relatively  high  data  rates,  and 
possibly  in  the  use  of  expensive  high  speed  memory. 
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S.3  COMMENTS  ON  FUTURE  RESEARCH 

This  investigation  has  raised  many  questions  which 
require  further  work  both  theoretically  and  experimentally. 

We  include  here  brief  descriptions  of  a  few  of  these  ques¬ 
tions  with  intuitive  comments  on  their  implications. 

5.3.1  Sample  Amplitude  Error 

The  velocity  samples  Ui  will  contain  "small"  errors 
from  sources  which  include  photon  and  electronic  noise, 
threshold  period  measurement  error,  ±  1  clock  count  error, 
and  inversion  round-off  error.  In  addition,  there  will  be 
occasional  "large"  errors  when  a  noise  burst  succeeds  in 
satisfying  the  check  logic  of  the  LV  processor.  Finally, 
there  is  a  different  type  of  effect  due  to  particle  lag, 
velocity  gradients,  and  finite  measurement  volume  size. 

As  a  first  approximation,  we  may  assume  that  the 
"small  errors"  collectively  form  a  zero-mean  wide-band  random 
noise  process,  n(t)  randomly  sampled  at  instants  t^.  If  n(t) 
is  independent  of  the  velocity  U(t)  then  the  actual  signal 

s(t)  =  n(t)  +  U(t) 

has  true  autocovariance  given  by 

C.Ct)  =  C(t)  +  C  «  CO 
s  n  ' 

where  CR  is  a  constant  and  600  is  the  unit  impulse  function. 
The  resulting  true  spectrum  is 

ssCO  *  SCf)  +  cn 
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i.e.,  the  value  of  Cn  is  the  white  noise  level  added  by  the* 
error  process  n(t) . 

A 

The  bias  error  effect  of  n(t)  in  the  estimate  C(k) 
will  be  confined  largely  to  C(0).  This  means  that  the  noise 
will  appear  in  the  measurement  process  as  turbulence  intensity 
which  is  not  actually  present.  We  observe,  however,  that  with 
the  simplified  model  assumed  for  the  noise  n(t)  the  value  of 
the  white  level  Cn  may  be  measured  independently  through  spec¬ 
trum  estimation  by  using  the  fact  that  the  true  spectrum  must 
go  to  zero  at  the  higher  frequencies.  Unfortunately  in  the 
widest  bandwidth  estimation  problems  where  the  ±  1  clock  errors 
are  most  likely  to  be  serious,  the  noise  spectrum  is  least 
likely  to  be  white.  The  immediate  result  will  probably  be  a 
small  signal- dependent  broad-band  spectral  bias  level  which 
will  add  further  difficulties  to  the  tail  estimation  prob¬ 
lem. 

The  second  class  of  errors,  ’’large  errors”,  occur 
relatively  seldomly,  but  the  much  larger  magnitude  implies 
that  these  errors  could  also  produce  a  non-negligible  addi¬ 
tive  white  spectrum.  Partial  remedies  include  the  data- 
check  circuits  and  preprocessing  steps  of  estimation  of  a 
and  throwing  out  excessively  large  or  small  estimates. 

These  practices  are  presently  being  incorporated  in  LV  ~ 
processors  and  will  aid  significantly  in  the  spectrum  esti¬ 
mation  problem. 

The  third  type  of  error  is  associated  with  the 
minor  particle  lag  and  velocity  gradient  effects.  These  are 
the  most  fundamental  in  that  they  will  be  present  even  if 
the  LV  system  is  error  free.  It  is  common  to  be  worried 
about  the  nontrivial  problem  of  the  low-pass  effect  of  par¬ 
ticle  inertia  on  spectrum  estimation.  We  have  left  that  one 

*  Further  analysis  has  revealed  that  this  white 
level  is  also  proportional  to  At.  See  "A  Discussion  of 
Limitations  and  Extensions  of  Power  Spectrum  Estimation 
with  Burst  Counter  LDV  Systems.",  W.  T.  Mayo,  Jr.,  to 
be  published  in  the  Proceedings  of  the  1974  Purdue  LDV 
Workshop . 
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for  the  particle  dynamicists .  A  more  subtle  effect. of 
minor  particle  lag  is  the  fact  that  two  particles  entering 
the  measurement  volume  at  different  velocities  and  closely 
spaced  in  time  will  appear  spectrally  as  part  of  a  wide¬ 
band  random  process.  The  same  type  of  effect  will  result 
from  velocity  gradients  when  particles  enter  opposite  sides 
of  the  probe  volume  closely  in  time.  These  effects  will 
not  only  increase  the  apparent  turbulence  intensity  but 
also  contribute  difficult- to-characterize  high-frequency 
spectral  error.  They  will  be  most  insidious  in  the  measure¬ 
ment  of  nearly  steady  flows  where  the  actual  turbulence 
levels  are  quite  low. 

A  related  effect  will  accompany  high  turbulence 
and  small  turbulence  scale:  while  the  finite  measurement 
volume  size  acts  as  an  averager  or  low-pass  filter,  it 
also  will  artificially  create  high  frequency  spectra  due 
to  random  timing  effects.  Thus  some  of  the  difficulties 
discussed  by  W.  George  [29]  under  the  assumptions  of  high 
density  seeding  and  Gaussian  statistics  also  exist  in  a 
different  form  in  the  single-particle  LV  systems. 

5.3.2  Cross-Correlations  and  Cross -Spectra 

We  wish  to  point  out  here  that  there  is  no  funda¬ 
mental  reason  for  not  extending  the  concepts  we  have 
developed  to  two-channel  cross-correlations  and  cross¬ 
spectra  and  that  this  is,  in  fact,  being  done  at  the  Lockheed- 
Georgia  Company,  In  this  regard,  however,  we  note  that 
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there  are  a  few  subtleties  in  measuring  cross-correlations 
between  two  separate  probe  locations  which  are  different 
than  simply  measuring  two  components  at  the  same  point;  and 
further  there  is  a  difference  in  a  two  component  measurement 
with  45°  crossed  fringes  and  a  one-component  measurement  with 
normal  fringes. 

First,  the  matter  is  straight-forward  for  two 
orthogonal  measurements  at  the  same  point  if  it  is  required 
that  both  channels  produce  "good  data"  simultaneously  and 
the  same  time  interval  measurement  applies  to  both.  The  only 
problem  would  seem  to  be  the  tolerances  on  the  "simultaneous" 
requirement.  However,  we  note  that  this  requirement  may  bias 
the  selection  of  "good  data"  toward  direction  vectors  along 
the  bisector  of  the  orthogonal  component  directions.  This 
results  when  Bragg- cell  systems  are  not  used  because  off- 
axis  direction  vectors  produce  a  higher  probability  of  "good 
data"  in  one  channel  than  in  the  other. 

If  simultaneity  is  not  required,  then  the  algorithms 
may  become  more  complex  as  they  will  in  separate-point  mea¬ 
surements  where  each  channel  has  a  partially  correlated  or 
independent  random  sampling  process.  The  correlation,  or 
lack  of  it,  between  the  sampling  processes  is  irrelevant  if 
they  are  jointly  independent  of  the  velocity  process.  The 
algorithms  for  time  difference  computations  are  more 
complicated  because  both  backward  and  forward  time  differ¬ 
ences  must  be  computed. 

Finally  we  note  that  with  space  correlations  involv¬ 
ing  zero  delay  in  time  the  concept  of  simultaneity  must  once 
more  be  considered.  Finite  tolerances  are  required  to  obtain 
any  data  points,  and  the  joint  data  rate  will  be  similar  to 
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that  of  one  element  of  the  autocovariance  histogram: 
generally  much  lower  than  the  individual  channel  data  rate . 

This  difficulty  cannot  be  overcome  by  indiscriminately 
opening  up  the  simultaneity  tolerance.  It  must  be  similar 
to  the  restrictions  on  At  in  the  autocovariance  estima¬ 
tion  to  prevent  the  averaging  to  zero  of  small  scale  compon¬ 
ents  . 

5.3.3  Correlations  Between  Mean  Sample  Rate  and  the  Velocity 
Process . 

In  informal  monthly  progress  reports,  we  derived 
the  fact  that  if  scatterers  have  a  uniform  Poisson  random 
distribution  in  volume,  then  the  expectation  A(t)  of  the 
sampling  process  is  not  constant  in  time.  Instead,  it  is 
proportional  to  velocity  and  the  estimate  of  the  mean 
produced  by  equation  (2.13)  is  biased  due  to  the  correla¬ 
tion  between  the  sampling  process  and  the  velocity  process. 

This  work  is  not  reproduced,  because  McClaughlin  and  Tieder- 
man  have  since  reported  similar  conclusions,  with  extension 
of  the  same  concepts  to  the  turbulence  intensity  estimate, 
in  detail  [30] .  We  note  that  subsequently  we  pointed  out 
that  the  bias  of  the  mean  estimate  vanishes  in  first  order 
approximation  when  the  average  of  the  burst  period  measure¬ 
ments  is  computed  before  inversion  and  scaling  to  velocity. 
Barnett  and  Bentley  [31]  have  extended  this  discussion 
which  has  now  become  a  topic  of  current  debate. 

It  is  not  unreasonable  to  suspect  that  correla¬ 
tion  between  velocity  fluctuations  and  mean  data  rate  will  also 


70 


AEDC-TR -74-53 


adversely  affect  the  spectrum  estimation  problem  under  high 
turbulence  intensity  conditions.  Unfortunately  we  are  now 
dealing  with  non- stationary  statistics  for  which  the  use  of 
interchangable  time  and  statistical  averages  is  not  appro¬ 
priate.  A  fallacy  underlying  the  whole  discussion,  however, 
is  the  lack  of  experimental  evidence  concerning  the  nature 
of  the  correlations  involved.  None  of  the  theory  mentioned 
thus  far  even  remotely  relates  to  the  interpretation  .of  jet 
mixing  data,  for  example,  in  which  the  mean  data  rate  depends 
on  which  air  source  each  small  volume  increment  originated  in. 
Considerably  more  definitive  work  is  needed  in  this  area. 
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SECTION  VI 

SUMMARY  AND  CONCLUSIONS 

The  objectives  of  this  research  have  been  accom¬ 
plished.  It  has  been  experimentally  demonstrated  that, 
with  randomly  timed  sampled  data,  power  spectra  may  be 
estimated  regardless  of  the  average  data  rate.  This 
means  that,  in  principle,  turbulence  power  spectrum  mea¬ 
surements  are  possible  with  unseeded  air  at  very  low  data 
rates  with  real-time  on  line  minicomputer  data  processing, 
without  the  use  of  large  high  speed  memory.  The  algorithms 
for  accomplishing  this  have  been  identified. 

The  principal  difficulty  lies  in  the  amount  of 
data  collection  time  required  to  determine  the  low  level 
spectral  tails  with  practical  error  tolerances.  Theory  has 
been  developed,  with  excellent  agreement  with  experiment, 
for  prediction  of  the  spectrum  estimate  error  in  terms  of 
the  algorithm  and  LV  system- flow  parameters.  The  results 
indicate  that  the  optimum  mean  data  rate,  for  efficient 
estimation  of  the  major  distribution  of  turbulence  power, 
is  twice  the  equivalent  power  bandwidth  of  the  spectrum; 

(not  twice  the  "highest"  frequency) .  However,  the  attainment 
of  practical  error  levels  in  spectrum  tail  estimation  may 
impose  very  stringent  requirements  including  low  resolution 
and  excessively  long  data  collection  and  processing  time. 

For  estimation  of  the  spectrum  tail,  the  most  significant 
error  and  time  reduction  is  accomplished  by  increase  of 
the  mean  data  rate  to  values  exceeding  twice  the  equivalent 
power  bandwidth  with  no  apparent  limit.  An  example  calculation 
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in  which  the  spectrum  tail  is  down  26  db  at  50  KHz  indicates 
that  a  mean  data  rate  of  23.5  K/sec  may  provide  practical 
measurement  times  and  error  levels.  In  the  same  example, 
real-time  continuous  minicomputer  processing  would  not  be 
feasible,  but  the  use  of  a  high-speed  disc  could  be  avoided 
through  recording  the  data  in  segments  with  gaps. 

A  preliminary  design  guideline  is  presented  and 
alternate  implementations  with  differing  tradeoffs  in  cost 
and  time  are  discussed.  At  the  present  time,  the  recommended 
implementation  would  use  a  minicomputer  as  a  buffer  and  pre¬ 
processor  for  recording  data  segments  on  digital  tape  or  for 
transmission  to  a  large  time- shared  computer.  Later,  when 
the  remaining  unknowns  which  have  been  identified  have  been 
more  fully  explored,  more  efficient  implementations  which 
destroy  the  data  as  it  is  obtained  may  be  appropriate. 
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APPENDIX  I 

TWO  SPECTRAL  ESTIMATORS  BASED 
ON  THE  PERIODOGRAM  APPROACH 

1.1  INTRODUCTION 

One  general  approach  to  the  estimation  of  broad¬ 
band  spectra  from  random  samples  is  to  approximate  the 
periodogram,  where  the  raw  periodogram  is  defined  as 

a  1  +t/ 2  2 

StO)=t  1/  u(t)  exp(-jwt)dt | 

1  1  -T/2  (I~l) 


with  u(t)  assumed  to  be  a  zero  mean  random  process  and  T 
is  the  record  length.  The  periodogram  estimate  is  made  to 
converge  by  averaging  many  estimates  from  short  segments  of 
data  [27] . 

1.2  ESTIMATOR  1 


The  most  straightforward  approach  is  to  use  a 
rectangular  approximation  of  the  raw  periodogram  integral, 
equation  (1.1).  After  some  manipulation,  the  result  is 


(1-2) 
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where 


Ai  "  ti  ti-l 


T  =  tN  ‘  t0 


1.3  ESTIMATOR  2 


Another  approach  is  to  first  use  the  periodo- 
gram  integral  to  estimate  the  spectrum  of  the  process 


Z2(t)  =  u(t)  Zx (t)  =  l  uA  6(t-ti) 

i=-<» 


+  OB 


(1-3) 


where  Z^t)  -J  aCt-t^, 


a=  -<*> 


It  is  assumed  Z^ft)  is  a  Poisson  impulse  process  with 
sample  rate  A.  With  this  estimate,  it  is  possible  to 
compute  an  estimate  of  S(w) ,  the  power  spectrum  of  U(t), 
because  it  is  known  that 


Sr(w)  *  A2  S(u)  +  Aa2  (1-4) 

where  Sr(w)  is  the  power  spectrum  of  Z2(t)  and  a2  is  the 
variance  of  U(t).  The  periodogram  estimate  of  Sr(u)  is 

ST(w)  =  f  |/  Z,(t)  exp  (- ju>t)dt  |  2 
1  o 

,  N  (I'5) 

■  ^\l  ui  expC-jwtj)  |2 

i=l 


78 


AED  C-TR-74-53 


and  the  estimate  of  S(o))  is 

*  1  f\  n  2 

s2Cu3)  -  J2  Sr^  "  — 


(1-6) 


For  real  computation,  it  is  necessary  to  use  estimates  of 
A  and  oz  which  are 


A  = 


A 


N 

T 


1  N  j 

-  n.E  ui 

i=l 


(1*7) 


Substitution  of  these  estimates  and  some  algebra  leads  to 


-  T  N  _  N 

S2  (w)  =  772  Kl  ui  cos  wti)Z+0!  u*  sin  wt.) 
N  i=l  i=l  1  1 


N 


(I~8) 


-I  «|] 

i  =  l 


1.4  EXPERIMENTAL  COMPARISON 


The  two  estimators  were  programmed  and  debugged 
using  both  uniformly  and  randomly  sampled  data  produced  by 
the  electronic  simulation  facility.  To  compare  the  esti¬ 
mators  when  used  with  the  same  data,  a  zero  mean  bandpass 
process  was  randomly  sampled  at  A  =  166  samples/sec.  The 
same  20  records  of  this  data  were  then  processed  with 
Estimators  1  and  2  and  with  the  autocovariance  estimator 
given  by  equation  2.33  of  the  text  (noted  as  Estimator  3 
in  Figure  1-3).  The  resulting  estimates  are  shown  in 
Figures  1-1,  1-2,  and  1-3  where  the  asterisks  represent 
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measurements  of  the  true  spectrum  as  made  with  a  wave 
analyzer, 

■ //.-Figure :  I- 4  is  the  result  of  reducing  M  by  a 
factor  of/4.;with  a  decrease  in  resolution  (which  was  four 
times  more  than  required  in  Figure  1-3)  to  Af  =  l/(2MAx) 
and  corresponding  reduction  of  estimate  variability. 

The' result  with  Estimator  1  was  not  satisfactory, 

-■  -  \’  -V 

while  the  results  with  Estimators  2  and  3  were  comparable. 
The  significant  differences  were  that  Estimator  2  required 
nearly  twice  as  much  computer  time  as  Estimator  3  and 
depended  heavily  on  the  stationarity  assumption  of  the 
sampling  process  for  the  prediction  of  the  white  level. 
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FHEBUENOf  -  HEKJZ 

Figure  1-1.  Spectrum  Estimate  by  Estimator  1 
(A  =  166  samples/sec) 

X  -  Estimate 

*  -  Scaled  Wave  Analyzer  Data 
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X  -  Estimate 

*  -  Scaled  Wave  Analyzer  Data 
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RB3LEH  -  HERT7 


Figure  1-3.  Spectrum  Estimate  by  Estimator  3 
(X  =  166  samples/sec) 

X  -  Estimate 

*  -  Scaled  Wave  Analyzer  Data 
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FREQUENCY  -  KERIZ 

-4.  Spectrum  Estimate  by  Estimator  3  with 

Additional  Smoothing ' (A-166  samples/sec) 
X  -  Estimate 

*  -  Scaled  Wave  Analyzer  Data 
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APPENDIX  II 


DERIVATIONS  OF  ESTIMATOR  ERRORS 


In  this  appendix  we  provide  the  details  concern¬ 
ing  the  random  sampling  bias  error  of  the  autocovariance 
function,  the  histogram  effects,  and  the  autocovariance 
and  spectrum  estimator  variability  errors. 


II. 1  RANDOM  SAMPLING  BIAS  ERROR 


Equation  (2.24).  for  C(0)  is  unbiased. 
(2.22)  for  C(k)  ,  k?*0,  may  be  expressed  as 


Equation 


i  H(k) 

C(k)  =  HtlT  £=1  Cm 


(II.l) 


where  C  is  a  u.u.  product  for  which 
m  i  j  r 

t  =  |t.-t.|  Cn-2) 

m  1  l  j  1 

satisfies  the  constraint  of  equation  (2.23).  Using  bracket 
notation  for  expectation  and  using  conditional  expectation, 
we  obtain 

(C(l°)  =  (bTO  Ji  Htk)  )Tin,  H(k)  UI-3) 

But  the  expectation  of  Cm  depends  only  on  the  velocity  process 

U(t)  and  the  time  difference  t  .  The  result  is 

m 

^irJrm’  =  ^CiJTm^  "  C^Tm^  fII-41 
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Taking  the  expectation  with  respect  to  rm  inside  the 
summation  gives 

V  y  i  H(k)  »  v 

<cm>  =<htst  £=1  LC(T)  pk(T)  dT>  hoo  cix^s) 

where  P^Ct)  is  the  conditional  probability  density  for  the 
value  of  t  .  given  that  t  lies  in  the  +At/2  interval  about 
kAt.  Since  the  integral  does  not  depend  on  H(k),  the  summa¬ 
tion  produces  H(k)  identical  terms  and  the  result  is  as 
given  in  (2.26) . 

The  function  P^CO  may  be  obtained  as 


Pk(T) 


Pn<t> 

kAx  +  Ax/2 
/  PN  (r)  dx 

kAr  -  At/2 


CII-6) 


where  p^Cx)  is  the  conditional  probability  density  function 

of  Tm  given  that  N  consecutive  samples  were  obtained.  The 

time  differences  t  may  be  divided  into  N-l  exclusive  sets 

according  to  the  order  of  the  time  difference:  there  are  N-l 

first  differences,  t.  -  t.  N-2  second  differences,  t.  -  t.  - 

1  l-l  l  i-2 

...;  and  one  N-l  difference,  t  -t, .  The  total  set  of  time 

n  l 

differences  is  the  set  union  of  N-l  mutually  exclusive  sets. 

We  may  therefore  express  the  density  p^Cx)  as  a  probability 
weighted  sum  over  the  respective  densities  for  first,  second, 
third,  etc.,  time  differences.  The  required  functions  f^(x) 
are  known  for  the  case  of  stationary  Poisson  sampling.*  The 
result  is 


PNW  ■  L  Pi  fi(T) 


OW] 


*Modifications  of  the  Erlang  functions.  See  page  558 
of  Papoulis  [21]. 
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where  =  probability  of  ith  order  difference  and 

\  r-  T^11 

f. 


C  X1  T1 

ijW  -  j  777777 


exp  [-Xt] 


t  >  0 
t  <  0 


cn-8) 


and  where  can  be  obtained  by  dividing  the  number  of 
ith  order  differences  (N-l)  by  the  total  number  of  dif 
ferences 


N-l 

I  (N-i) 
i-1 


N2-N 

T~ 


CII-9] 


with  the  result 


2  (N-  i 

i  "  N2.N  CII-10) 


Finally  we  see  that  an  exact  expression  for 
pk(x)  may  be  obtained  as  a  finite  sum  of  integrals  in¬ 
volving  the  f i C T 3  functions  if  (11-10),  (II-8),  and  (II-7) 
are  substituted  in  (II-6) .  The  result  is  a  complicated 
expression  which  could  be  approximately  evaluated  by  power 
series  or  numerical  integration; . but ■ this -is  not  yet  been 
undertaken. 


II . 2  Histogram  HCk) 


Using  the  same  type  of  derivation  as  in  the  pre¬ 
vious  section  we  would  obtain 


H(k) 


N2-N 

~2 


kAx  +  At/2 
/  pn(t)  dx 
kAx  -  Ax/2 


CII-II] 
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This  expression  has  not  been  evaluated  and  compared  with 
experiment. 

Using  the  approximate  alternate  formulation  as 
in  equations  (2.25)  with  the  conditions  for  (2.7),  we  ob¬ 
tain  the  result  given  in  equation  (2.27)  after  noting  that 

Iitk>  (N'-k)  p2  CII-1Z) 

by  replacing  N*  =  T/At  =  N/(XAt).  The  extent  of  validity 
of  the  approximate  result  is  not  known,  but  good  agreement 
with  experiment  was  obtained  in  the  experiments  performed. 


II. 3  VARIANCE  OF  THE  AUTO COVARIANCE  ESTIMATE 


In  addition  to  the  assumptions  of  Section  2.1.1 
several  more  assumptions  are  used. 

a)  At  is  small  enough  to  prevent  bias  error 

b)  RMS  deviation  of  H(k)  is  small  (large  N) 

c)  H(k)  =  constant,  k  <  M  (large  N) 

d)  individual  product  pairs  are  independent  events 

A 

e)  the  true  mean  is  known  (U  =  U) 

Assumption  d)  is  not  strictly  ever  valid;  however,  for  suf¬ 
ficiently  low  values  of  X  there  are  negligibly  few  instances 
in  which  two  product  pairs  with  the  same  lag  value  occur  close 
enough  together  in  time  to  be  appreciably  correlated.  It  was 
initially  assumed  that  X/2B  <<  1  would  be  required  to  obtain 
d) .  The  question  has  not  been  resolved  analytically,  but  the 
following  results  appear  to  be  valid  experimentally  for  values 
of  X  up  to  2B  (and  even  higher  for  the  stopband  frequency  com¬ 
ponents  . ) 
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A 

Under  our  assumptions  C(k)  is  unbiased  so  that 


<feT(k)>  =  <"cCk)  -  C(k)>  =  0 


Cl  1-13) 


and  the  mean  square  error  is 

v  ,  H(k)  H(k)v  Vv 

ackW  1<CU  CJV\S>  -  <=*00 

m-i  n-i  ^H(k) 


CII-14) 


with  and  each  being  product  pairs,  as  in  (II-l),  which 
are  independent  of  H(k) .  Expanding  one  term  of  the  summation 
we  have  H(k)  terms  where  m  =  n: 

<V(t)  u2(t  +  Tm)  It^,-  a4  +  2  C2(xm)  CII-15] 


2 

and  H  (k)  -  H(k)  terms  m  i  n  which  with  the  assumed  indepen¬ 
dence  gives 


^u(t)  u(t+Tm)  u(t')  u(t,+xn)^ 

"  ^f(t)  UCt+Tni),>  U(t’+Tn)^> 


C(t  ) 
v  nr 


C(t  ) 
n' 


(11-16) 


Assumption  a)  allows  replacement  of  xm,  by  kAx 
with  negligible  effect  so  that  (11-14)  becomes 


ack=(iFW  {H(k)  [°"  +  2c2(k)l  + 


[H2 (k)  -  H (k) ]  C2 (k) }  )H(k)  -  C2  (k) 


Cir-i7) 
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We  observe  that  equation  (2.31)  results  if  the  RMS  deviation 
of  H(k)  is  small  so  that 


1 


(11-18) 


II. 4  VARIANCE  OF  THE  SPECTRAL  ESTIMATE 


With  the  additional  assumption  that  MAr  is  large 
enough  so  that  negligible  smoothing  bias  error  exists,  the 
error  e^(i)  in  the  spectral  estimate  at  iAf  is  zero  mean. 

The  mean  square  error  is  obtained  by  inserting  equation  (11.13) 
in  equation  (11-33)  squaring  the  result,  expanding  as  a  double 
sum  as 


a£(i) 


*  4At: 


M 

1 

=■0 


where 


M 

I  ' 

n=0 


w 


(m)  w’ (n)  (eT(m)  eT  (n)  )  cos  (H-19) 


cos 


r— i 

1  M  J 


w' (k)  =  w(k),  k/0;  w' (k)  =  1/2,  k=0 


(11-20) 


But  since  S(I)  is  unbiased  and  eT(m)eT(n)  =0  for  m^n 
because  of  the  assumed  independence  the  result  becomes 

o|(i)  =  2At2  l  w’2(k)  a2^k)  (1+  cos  ) 

k=0  (11-21) 
2  2 

neglecting  the  fact  that  w,4(0)  =  1/4  and  that  w  (0)  =  1, 
equation  (2.41)  results  from  inserting  equation  (2.27)  into 
equation  (2.31),  with  (2.28)  satisfied,  then  using  equation 
(2.31)  in  (11.21) . 
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APPENDIX  III 

EXPERIMENTAL  FACILITIES  AND  DATA  PROCESSING 

11 1. 1  ELECTRONIC  SIMULATION  FACILITY 

111. 1.1  Introduction 


The  electronic  simulation  facility  used  the 
random  emissions  from  a  radioactive  source  to  control 
an  analog-to-digtal  converter  (ADC)  and  a  special  purpose 
time -interval  counter.  The  ADC  produced  sample  magnitudes 
and  the  time- interval  counter  counted  the  number  of  clock 
pulses  between  samples  to  indicate  the  time  between  conse¬ 
cutive  samples.  The  digital  data  from  the  ADC  and  the 
time- interval  counter  was  read  by  a  minicomputer  and  stored 
on  magnetic  tape  for  off-line  processing.  The  system  block 
diagram  in  Figure  III-l  shows  the  general  arrangement  of 
equipment.  Each  component  is  discussed  below. 

III. 1.2  Radioactive  Source  and  Proportional  Counter 

The  radioactive  source  was  a  laboratory  sample 
of  carbon  14  of  .001  yc  activity.  The  gas  proportional 
counter  detected  the  emission  of  alpha  and  beta  particles 
and  generated  a  -1.5  v  pulse  with  duration  1  ys  for  each 
detection  event.  The  dead  time  between  events  was  of  the 
order  of  a  few  microseconds.  The  number  of  events  per  unit 
time  was  assumed  to  have  a  Poisson  distribution  with  mean 
rate  A.  The  mean  rate  was  controlled  by  covering  the 
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sample  with  various  thicknesses  of  aluminum  foil,  while 
maintaining  a  constant  voltage  on  the  proportional  counter 
electrode  of  approximately  2000  volts.  The  mean  rate  was 
stable  under  these  conditions. 

111. 1.3  Control  Logic  and  Counter 

This  component  contained  three  sections: 

1)  circuitry  to  interface  the  output  pulse  of 
the  gas  proportional  counter  to  TTL  digital  logic, 

2)  TTL  logic  to  generate  control  pulses  for  the 
ADC  and  the  time- interval  counter, 

3)  a  time- interval  counter  consisting  of  standard 
TTL  packages  and  a  Schmitt  trigger  clock  to  drive  the  counter. 

The  time- interval  counter  was  16  bits  (15  magni¬ 
tude  bits  and  a  sign  bit).  The  sign  bit  was  included  solely 
for  convenience  in  interfacing  with  Hewlett-Packard  logic. 

The  time  resolution  of  the  counter  was  approximately  2  us 
(the  period  of  the  clock)  and  the  upper  range  was  (21S-1)  x 
2  ys  =  .066  seconds. 

111. 1.4  Noise  Generator  and  Filter 


The  noise  generator  provided  a  zero  mean  Gaussian 
signal  which  was  approximately  "white"  over  the  range  5  Hz 
to  20  KHz.  The  filter  was  an  adjustable  fourth  order 
Butterworth  active  filter  that  could  be  used  in  a  lowpass, 
bandpass  or  highpass  mode.  The  3  db  cut-off  frequencies 
were  variable  from  10  Hz  to  1.1  MHz.  The  roll-off  rate  was 
24  db/octave.  In  the  lowpass  mode,  experimental  measurements 
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with  a  wave  analyzer  showed  that  the  response  of  the 
filter  was  not  flat  below  approximately  20  Hz.  This  was 
verified  in  the  spectral  estimates.  The  wave  analyzer 
measurements  also  showed  that  the  cutoff  frequencies  indi¬ 
cated  by  the  filter  dial  settings  were  low  by  5%.  Figure 
III- 2  shows  a  plot  of  wave  analyzer  data  (dashed  line)  and 
a  plot  of  the  theoretical  Butterworth  curve  (solid  line) 
for  comparison. 

I I I. 1.5  DC  Source  and  Summer 


The  DC  Source  was  simply  a  DC  power  supply.  The 
DC  level  was  added  to  the  output  of  the  filter  to  provide 
a  non- zero  mean  signal.  To  perform  the  addition,  a  summer 
was  built  with  a  general-purpose  operational  amplifier. 

The  summer  provided  a  gain  of  approximately  100.  Adjust¬ 
ment  of  the  DC  power  supply  and  the  output  level  of  the 
noise  generator  provided  a  nonzero  mean  signal,  with  variable 
turbulence  intensity  and  known  spectrum,  in  the  input  range 
of  the  ADC. 

III. 1.6  Analog  to  Digital  Converter 

The  ADC  had  an  input  range  of  ±  5  volts  and  the 
digital  output  consisted  of  nine  magnitude  bits  and  a  sign 
bit  representing  a  decimal  range  of  -511  to  +512.  Negative 
magnitudes  were  in  two's  complement  form.  The  ADC  required 
a  12  ps  trigger  pulse  which  was  generated  by  the  control  logic 
described  earlier. 

The  minimum  sampling  interval  was  approximately 
60  ys,  consisting  of  50  ys  dead  time  for  the  ADC  (maximum) 
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and  a  10  ys  allowance  for  execution  time  of  computer 
instructions.  Randomly  spaced  samples  occurring  at  smaller 
time  intervals  were  ignored. 

III. 1.7  Minicomputer 


The  data  was  collected  with  a  Hewlett-Packard 
2114B  minicomputer  (8K  memory)  and  recorded  on  a  Hewlett- 
Packard  7970B  Digital  Tape  Unit  in  integer  format  in  data 
segments  consisting  of  500  words  (250  samples)  per  segment. 
An  existing  program  (ADCMT)  written  at  Texas'  A§M  was  modi¬ 
fied  for  use  with  the  ADC  and  the  time  interval  counter. 

III. 1.8  Instrumentation 

A  laboratory  time- interval  counter  was  used  to 
count  the  number  of  random  pulses  in  a  given  interval  to 
determine  the  average  sample  rate.  A  laboratory  wave 
analyzer  was  used  to  measure  the  distribution  of  power  vs 
frequency  at  the  output  of  the  analog  filter.  The  magni¬ 
tude  of  the  time-varying  component  was  measured  with  a 
true  pms  voltmeter. 

III. 2  DATA  PROCESSING 

III. 2.1  Data  Files 


The  data  for  each  experimental  run,  recorded  on 
a  Hewlett-Packard  Minicomputer  system,  was  taken  to  an 
IBM  360/65  facility  (with  IBM  3420  tape  drives)  for  transfer 
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to  a  Master  Data  Tape.  The  systems  were  compatible  for 
integer  data,  except  for  word  length  (16  bits/word  for 
Hewlett-Packard;  32  bits/word  for  IBM  360/65).  Each 
experimental  run  was  assigned  a  separate  file  on  the 
Master  Data  Tape. 

III. 2. 2  Software 


Programs  were  written  in  Fortran  IV  to  batch 
process  the  experimental  data  in  data  segments  of  500 
words/segment,  where  1  word  =  16  bits.  The  autocovariance 
estimator  described  in  2.3.1,  eq.  (2.22),  was  implemented 
as  a  subroutine  by  keeping  running  accumulations  in  the 
arrays  H(k)  and  SUM(k)  until  a  predetermined  number  of 

A 

data  segments  had  been  processed.  The  estimate  C(k)  was 
then  formed,  followed  by  the  spectral  estimate  S(iAf) 
computed  from  eq.  (2.33).  The  spectral  estimation  sub¬ 
routine  used  is  reproduced  on  the  following  pages. 

Other  subroutines  were  written  to  read  data 
segments  from  the  Master  Data  Tape,  to  compute  the  -sample 
time  of  each  sample  assuming  t  =  0  at  the  beginning  of 
each  data  segment,  to  compute  segment  mean,  and  to  plot 
the  various  arrays  and  estimates  on  a  line  printer. 
Execution  times  of  the  order  to  45  sec  were  encountered 
when  processing  20  data  segments  with  this  package,  but 
this  could  be  decreased  considerably  by  reprogramming 
with  use  of  integer  arithmetic  operations. 
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III. 2. 3.  Typical  Output 

Figure  III-3,  a  and  b,  illustrate  the  difference 
between  a  linear  plot  and  a 'logarithmic  plot.  Because  of 
the  symmetrical  type  of  error  produced  by  the  estimator 
of  Section  2.4,  (not  positive  definite)  linear  plot  routines 
have  been  used  during  this  study.  The  two  figures  illus¬ 
trate  the  difference  in  appearance  for  normalized  RMS  error 
level  of  -17  db. 
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1 1 1 - 3  a.  Linear 
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III-3  b.  Logarithmic 


Figure  1 1 1 - 3 .  Linear  and  Logarithmic  Plots  of  a  Spectrum 
Estimate . 
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S UB ROUT IME  R ASPEC  (NUM, NUM 1 »M»  DEL TAU.AMPL tTtSUM»H,C ,PS, NF  » 


1KLM.KGM, JMAX, STEP! 

MODEL  5  REENTRANT  MOO  5  JANUARY  28» 1974 


SPECTRUM  ESTIMATION 

FROM  RANDOMLY  SAMPLED  RANDOM  PROCESS 

. NUM  -  NUMBER  OF  DATA  WORDS  ' iN"EACH  'RECORD 

NUM1  -  NUM  +  1 


M  -  NUMBER  OF  DELTA  TAU  SLOTS  USEO  IN  HISTOGRAM 
DELTAU  -DELTA  TAU  -  TIME  LENGTH  OF  EACH  SLOT 
"  AHPL  -  INPUT  DATA* AMPLITUDE  MATRIX 

(RECORD  MEAN  ALREADY  SUBTRACTED  __  _ 

T IME  MATRIX  -  ELAPSED  TIME  CORRESPONDING  TO  EACH  AMPLC I ) 
MATRIX  OF  SUBTOTALS  OF  AMPL ITUOES  CORRESPONDING  TO  EACH 
HISTOGRAM  SLOT 

MATRI X  OF  SUBTOTAL  NUMBER  OF  ELEMENTS  IN  EACH  SLOT 

CORRESPONDS  TO  SUM  MATRIX . 

AUTOCOVARIANCE  MATRIX _ 

POWER  SP  ECTRUM  MATR I X"  OF  ,;L  ENGT 

NUMBER  OF  FREQUENCIES  T,o"  BE  EVACUATED  FOR  -POME  ft  SPECTRUM 


KLM  -  NUMBER  OF  LAGS  LESS  THAN  OF  EQUAL  TO  M*>~ 

KGM  -  NUMBER  OF  LAGS  GREAfER  TrHAN’.M 

C  ONf  I NUE  .  . . "Cs . 

THE_FULL  SUBROUTINE  MUST  BE  X  ALLE  D4F  I RST  TO .  INITIAL  I ZEALL  STORAGE 
LOCATIONS.  "  CONTROL"  IS’ ''RETURNED'  "IMMEDIATELY  f  •  v  - 

PSTEP2  IS  CALLED  EACH  TIME^FOR  SOME. NUMBER 'OF  RECORDS, OFA INPUT 


CONTROL  IS  RETURNED  AFTER  THE  HISTOGRAM  AND  SUM/NATRI  X,-‘ARE„  .COM¬ 
PLETED  FOR  THAT  RECORD  '*  ’  j r  W'.„.  , 

P STEP 3  I  S  CALLED  TO  CREATE-  THEr  d MATRIX  FROM  THE  ACCUMULATED  SUM 
AND  H  MATRICES.  RETURNS  C*  f.  „  1  ? 


PSTEP4  MULTIPLIES.  C  BY  THE  W  RAMP.  TTHE  W.  RAMP  ^STARTS  AT*  l  AND 
FALLS  TO  0  IN  IH  STEPS  AND  REMAINS,  0  FOR  THE-'  RES f  OF  TH/C  MATR IJC 
PSTEP5  CREATES  THE  POWER  SPECTRUM  ESTIMATE .FOR. NF  FREQUENCIES 

I  START  -  STARTING  FREQUENCY  DIVIDED' BY  2.5 _ _  "" 

STEP  -  MULTIPLE  OF  2.5  TO  BE  STEP -IN  FREQUENCY,  . 

POWER  SPECTRUM  IS  RETURNED  AND  SUBROUTINE  'IS  ' TERMINATED 


uu.o  ;  I  .  !  uuu 


INTEGER  JMAX.HIM) ,C0, STEP 
_ PEAL  SUM(Hl,C<H),PStNF).TINUHl) .AMPL(NUM) 

...  jMTiAlIZE  CONSTANTS  AND  MATRICES  ' 

SUM0=0. 

H0=0 
KLM=0 
KGM=0 

DO  5  1=1, M 
"SUM(T)=o: 

5  H( 11=0 

return . . 

_ ENTRY  PStEP2  (NUN,AHPL,T, SUM,H,W,NUH1,KLH,KGH, JHAX) 

PRODUCE  SUM  AND  H  ARRAYS,  KEEPING  TRACK  OF  KLM  AND  KGM 

O  CO  10  1=1, NUM 

M  "10”  S  UV0="SU  MO+AMPL  (  1 1 +AMPL  ( I  } 

HO=NUM+HO 
DO  20  J=l , JMAX 
N  P=NUM— J 
DO  30  I=1,NP 
11=1+1 

DELTA =TC  I  i+J  ) -T (  ID" . 

K=DELTA* 1 1 ./DELT  AUl +1 • 5 
!F(K,GT.M)GO  TO  25 

29  H  ( K.)  =H(  K )  +1 

- STTMT  K)  =^UM(Kf+AMTsUfT+  3  J+AMPLin 

KLM=KLM+1 

- GTTTO'TO'  ■ 

25  KGM=KGM+1 

30  CONTINUE 
20  CONTINUE 

SON!  it^sumu 
HC 1 1  =  H0 
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oooo  \  o  o  o  «  o  o  o  o  o  r> 


RETURN 

ENTRY  PSTEP3  <SUM*H.C,M)_ 
PRODUCE  AUTOCOVARIANCE  MATRIX 


MM1=M— 1 

DO  40  K=2  *MM1 

IFlHtKJ .NE.01  GO  TO  40 

IF  HI  K)  =0  ♦  THEN  STRAIGHT  LINE  INTERPOLATE  BETWEEN  HIK-1I  AND  HIK+1) 

SUM(K)=SUM(K— U+SUMIK+ll 
H(K)=H(K-1  )  4-H  I  K+l ) 

40  CONTINUE  _ 

'"'TFT H(  M)  «’NE  #6T  GO  TO'  46 

IF  HlM)=Ot  THEN  HIM)=H(M— 1) 

SUM(M»=SUMIM-1) 

HIM} =H( M-l) 

"46"0U~WT=1TM . 

50  C(K|=SUM(K)/H(K| 

RETURN  •' 

ENTRY  PSTEP4  (CvIHfNI 

°R00UCE  JW  RAMP  IBARTLETT  WINDOW)  AND  MULTIPLY  C  MATRIX  BY. THIS 
WEIGHTING  FUNCTION 

00  60  K=I » IH 

W=I. -FLOAT  (K-D/FLOATUH) 

60  C(K)=W*CIK) 

I  HI  =  IH  +  I 

DO  61  k=IHltM 

61  C (K ) =0. 

RETURN 

ENTRY  PSTEP5  I NF , ISTART, STEP. PS) 
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E.S I!M ATE  POWER  SPECTRUM  FOR  DESIRED  FREQUENCIES 

DO  TO  J=I,NF 
I  =1 START+STEP*(J-1) 

SUMT  =  0« 

DO  80  K=2  » M 

S  UM  T  =  S  UMT +CJ  K )  *C  OS  (_3  .  1 4 1 5  9*  I  *  F  L  C  A  T  ( K- 1)  /  F  L  0  AT  |  M  )  ) 
SUMT=SUMt*SUMf+CI I)  ' 

PS( J)=DELTAU*SUMT 
RETURN 


